EXISTENCE AND STABILITY OF VISCOUS SHOCK PROFILES FOR 2-D 
ISENTROPIC MHD WITH INFINITE ELECTRICAL RESISTIVITY 



Abstract. For the two-dimensional Navier-Stokes equations of isentropic magnetohydrodynamics 
(MHD) with 7-law gas equation of state, 7 > 1, and infinite electrical resistivity, we carry out a 
global analysis categorizing all possible viscous shock profiles. Precisely, we show that the phase 
portrait of the traveling-wave ODE generically consists of either two rest points connected by a 
viscous Lax profile, or else four rest points, two saddles and two nodes. In the latter configuration, 
which rest points are connected by profiles depends on the ratio of viscosities, and can involve Lax, 
overcompressive, or undercompressive shock profiles. For the monatomic and diatomic cases 7 = 5/3 
and 7 = 7/5, with standard viscosity ratio for a nonmagnetic gas, we find numerically that the the 
nodes are connected by a family of overcompressive profiles bounded by Lax profiles connecting 
saddles to nodes, with no undercompressive shocks occurring. We carry out a systematic numerical 
Evans function analysis indicating that all of these two-dimensional shock profiles are linearly and 
nonlinearly stable, both with respect to two- and three-dimensional perturbations. For the same 
gas constants, but different viscosity ratios, we investigate also cases for which undercompressive 
shocks appear; these are seen numerically to be stable as well. 
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1. Introduction 

In this paper, we continue the investigations of [SZl EH IMaZ3l IMaZ4l iRal iRZl IZ5l ITZl IBiSZl 



[BFTT lBr2l IBrZl IHuZ2l IBHKZl iHLZl |HLyZl| |HLyZ2| IBHZ] on stability and dynamics of large 



amplitude viscous shock profiles, examining classical Lax-type and nonclassical overcompressive 
and undercompressive shocks occurring in isentropic magnetohydrodynamics (MHD) with infinite 
electrical resistivity. 

Existence of large-amplitude profiles for full (nonisentropic) magnetodydrodynamics was studied 
in pioneering works of Germain and Conley-Smoller (21 ICS1[ ICS2| . making use of properties of 
the traveling- wave ODE as a gradient system and of Conley index techniques. Further investiga- 
tions have been carried out by Freistuhler-Szmolyan |FSj using geometric singular perturbation 
techniques and by Freistiihler-Roehde |FR1| IFR2j using a combination of bifurcation analysis and 
numerical approximation. In this generality, the traveling-wave ODE for MHD profiles is a six- 
variable dynamical system, with up to four rest points corresponding to endstates of various inviscid 
shock waves. For an ideal gas law, it is known that fast and slow Lax shocks always possess a vis- 
cous profile. In certain special cases, or in certain limiting ratios of viscosity, heat conduction, 
etc., it is known that intermediate shocks do or do not possess profiles; however, in general, the 
profile existence problem for the full nonisentropic case is accessible at present only numerically. 
For further discussion, see |FS1 IFR1| IFR2| and references therein. 

In the present work, we examine in detail the restricted case of isentropic flow with infinite elec- 
trical resistivity, in two dimensions, for which the traveling-wave ODE becomes a planar dynamical 
system. This example exhibits the main features of the general case, in a simpler setting conducive 
to systematic numerical investigation. 

Specifically, for a rather general equation of state (convex, decreasing in specific volume, and 



blowing up at least linearly with density as density goes to infinity) we show in Sections 2.43 and 



(4.1 ) that the phase portrait of the traveling- wave ODE generically consists of either two rest points 
connected by a viscous Lax profile, or else four rest points, two saddles and two nodes. In the latter, 
four rest point configuration, the Lax shocks involving consecutive rest points ordered by specific 
volume always have connecting profiles. The remaining, "intermediate" shocks may or may not 
admit profiles, depending on the ratio of parallel to transverse viscosity. Specifically, we show in 



Section 4.3 by phase plane (and, separately, by singular perturbation) analysis that, similarly as 
in the nonisentropic case |FS[ IG] , any intermediate shock with decreasing specific volume permits 
a connection for some viscosity ratios and not for others. By entropy considerations, shocks with 
increasing specific volume never have connecting profiles. Here, and elsewhere, we without loss 
of generality restrict discussion to the case of a left-moving shock. (For right-going shocks, the 
ordering would be reversed.) 

We supplement this abstract existence discussion by a systematic numerical existence study for 
specific parameter values in physical range. For the most common cases of monatomic or diatomic 



gas, 7 = 5/3 or 7 = 7/5, with standard viscosity ratio for a nonmagnetic gas (see (2.3 )), we find that 
there occurs only one profile configuration, with the nodes connected by a family of overcompressive 
profiles (intermediate shocks) bounded by Lax profiles connecting saddles to nodes in a four-sided 
configuration (one pair of opposing sides corresponding to slow and fast Lax connections, the other 
to intermediate Lax connections). Undercompressive profiles do not seem to occur in this parameter 
range. 

Next, restricting to the same parameters 7 = 5/3,7/5, and standard viscosity ratio, we carry 
out numerically a systematic stability analysis of these waves, using the general numerical Evans 



function techniques developed in [BUI iBrZl IHuZ2l iHLZl |HLyZl| |HLyZ2| lBHZllZ5] . Our results 



carried out up to extremely high Mach number (typically Mach 20 — 40, but in some cases up 
to Mach 10,000), indicate that all of the above profiles, both Lax- and overcompressive type, 
are spectrally stable in the generalized Evans function sense defined in jZHt IMaZ3j . both with 
respect to two-dimensional and three-dimensional perturbations. These results are described in 
Sections |5] and [7| By the abstract framework established in (MaZ3, M aZ4| [ZT] \R&\ IRZ] . this implies 
linearized and nonlinear time-asymptotic orbital stability, as described for completeness in Section 



2.5 Varying the viscosity ratio, we carry out case studies also for examples of undercompressive 
profiles. Numerically, these are seen to be (Evans, hence linearly and nonlinearly) stable as well. 
Finally, in Section [8] we discuss our results and suggest directions for further study. 

2. Preliminaries 

2.1. Equations and assumptions. In Lagrangian coordinates, the equations for compressible 
isentropic magnetohydrodynamics (MHD) take the form 

vt — u\ x = 0, 

uit + {p + (1/2W))(£ 2 2 + Bj)) x = (((2/i + rj)/v)u lx ) x , 
u 2t - {{1/hq)IB 2 ) x = {{ll/v)u 2x ) x , 

U3t ~ ((l/Mo)^)^ = {{^/v)u 3x ) Xl 

(vB 2 ) t - (Iu 2 ) x = ((l/o-fj, v)B 2x ) x , 
(vB 3 ) t - (Iu 3 ) x = ((l/afj, v)B 3x ) x , 



(2.1) 



where v denotes specific volume, u = (u\, u 2 ,u 3 ) velocity, p = p{v) pressure, B = (I,B 2 ,B 3 ) 
magnetic induction, / constant, and \x > and 77 > the two coefficients of viscosity, /j,q > the 
magnetic permeability, and a > the electrical resistivity; see |7Q £3 [J] for further discussion. 

We restrict mainly to the case of an ideal polytropic gas, in which case the pressure function 
takes form 



(2.2) 



p{y) = av 
3 



where a > and 7 > 1 are constants that characterize the gas, the limiting case 7 = 1 correspond- 
ing to the barotropic, or constant-temperature approximation and 7 > 1 corresponding to the 
isentropic, or constant-entropy approximation, of the ideal pressure law p(v, e) = Tv~ 1 e. Though 
we do not specify r], we have in mind mainly the ratio 



(2.3) 



77 = -2/z/3 



typically prescribed for (nonmagnetic) gas dynamics |Baj . (By rescaling space and time, we can 
rescale all transport coefficients by a common factor; the ratio rj/fi, however, is invariant.) 

In the thermodynamical rarified gas approximation, 7 > 1 is the average over constituent particles 
of 7 = (JV + 2)/N, where N is the number of internal degrees of freedom of an individual particle, 
or, for molecules with "tree" (as opposed to ring, or other more complicated) structure, 



(2.4) 



7 



2n + 3 
2n + P 



where n is the number of constituent atoms [Ba]: 7 = 5/3 ~ 1.66 for monatomic, 7 = 7/5 = 1.4 
for diatomic gas. 

An interesting subcase is the limit of infinite electrical resistivity a = 0, in which the last two 



equations of (2.1) are replaced by 

(2.5) (vB 2 )t - (Iu 2 ) x = 0, (vB 3 )t - (In 3 ) 3 



0., 



and only the velocity variables u = (ui,u 2 ,Us) experience parabolic smoothing, through viscosity. 
We can restrict further to the two-dimensional case, setting 113, B% = and dropping these variables 
from consideration, as we shall do for most of our investigations. 



2.1.1. Eigenvalues of the 2-d inviscid system. The inviscid version of system (2.1) in dimension 
two, is, introducing the scalar quantities B = B 2 , w = u 2 , 



(2.6) 



or, in quasilinear form, 



0, 

B 



u t +p x + j^B x 
w t - ±B X = 0, 



mo 

B t + f u x 



v J 



(2.7) 



where 



(2.8) 
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p'{v). This system has four eigenvalues of the form ±y / r±, where r± are the roots of 



c 2 + 



P + B\ I 2 . 
r+ c 2 



0. 



As the discriminant of (2.8) is positive for the two roots r + and r_ are positive real, verifying 

hyperbolicity. When B = 0, the discriminant can be zero for c 2 = We do not explicitly require 
this computation in our analysis, but include it for general interest /orientation. 
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2.2. Viscous shock profiles and the rescaled equations. A viscous shock profile of (2.1) is a 
traveling- wave solution, 

(2.9) (v, u, B){x, t) = (v, u, B)(x - st), 

moving with speed s and connecting constant states 



(2.10) 



(v±,u±,B±)= lim (v,u, B) (z) 

z— >±oo 



Such a solution is a stationary solution of the system of PDEs 
(2.11) 

where we have denoted u := (u 2 , 1*3), B := (B 2 , -B3), i.e., a solution of the system of ODEs 
(2.12) 



v t - sv x - u\ x = 0, 

uit ~ su lx + (p + (l/2fj, )\B\ 2 )x = (((2/U + r])/v)u lx ) x , 
Ut - su x - {{l/po)IB) x = ((fj,/v)u x ) x , 
k (vB)t - s{vB) x - (Iu) x = ((l/o-(j, v)B x ) x , 



' -sv' - u[ = 0, 
- S < + (p + (1/2/io)^! 2 )' = (((2/x + v)/vKY, 
-su> - ((1/im>)IB)' = ((»/v)u>)>, 

I _ s{v b)> - (my = {{i/a^v)&y. 



Integrating, we obtain 
(2.13) < 



— sv — u± = Ci, 

-sui + {p + {l/2fi )\B\ 2 ) = (((2// + ri)/v)u[) + C 2 , 
-su - ((l//io)/#) = ((/*/«)«') + C 3 , 
I - s (vB) - (Iu) = ((l/an v)B') + C 4 



for some constants of integration C := (C\, . . . , C4). 

For fixed C, the rest points of (2.13) comprise the possible endstates (v± , u± , B± ) that can 
be connected by a viscous profile with speed s, which necessarily satisfy the Rankine-Hugoniot 
conditions 



-s V = u 



(2.14) 



-s[ui] 



-s u 



p + 
B 

t'Q 



2/U 



-s[vB] = I[u] 

determining pairs of states connected by an inviscid shock wave, where 

[h] := h(v + ,u + , B + ) — h(v-,U-,BJ) 
denotes jump in the quantity h across the shock. 

2.2.1. Rescaled evolution equations. Following [HLZ { [HLyZl HLyZ2| [BHZ | , we now rescale 

'v u , o B ae _7_1 ' 



(2.15) 



(v,u,fiQ,x,t,B,a) 



es 



, e/io, —es(x — st), es 2 t, 



holding /t, a fixed, where e := V—, transforming (2.1), (2.2) to the form 



(2.16) 



v t + v x - u lx = 



1 

2/iQ 



IB 2 
IB 3 



U2t + U 2x 
U3t + U3x 

(vB 2 ) t + (vB 2 ) x - (Iu 2 ) 

(vB 3 ) t + (vB 3 ) x - (Iu 3 ) 
where p(v) = aw~ 7 . There is no change in [i or r\. 



(2M + ,)p) 

V V lx 

(u 2x \ 

V V ) x 



/' 
/' 



fu 3x \ 

V V I x 
\<J\XqV 

\au v 



B 2x 
B 3x 



2.2.2. Rescaled profile equations. Viscous shock profiles of (2.16 ) must satisfy the system of ordinary 
differential equations 

r v' - u\ = 0, 



(2.17) 



u' x + (p + (1/2^)\B\ 2 )' = (((2 M + 77)^X1)', 
v! - {{l/^)IB) ! = ((ji/vW, 
I (vB)' - (Iu)' = ((l/a^v)B')', 

together with the boundary conditions 

(v,Ui,u,B)(±oo) = (v,u 1 ,u 2 ,B)±. 

Evidently, we can integrate each of the differential equations from — oo to x, and using the boundary 
conditions (in particular V- = 1 and u_ =0), we find, after some elementary manipulations, the 
profile equations (after having introduced shorthand notation u = u±, w := u, B = B): 



(2.18) 
(2.19) 
(2.20) 

with u = v — 1. 



{2fi + rj)v' = v(v - 1) + v(p - p_) + — (B 2 - B 2 ), 

2^o 

/ y I 
fiw = vw (B — -B-), 

Mo 

B' = v 2 B — vB_ — Ivw, 



2.2.3. The case a = oo. When a = oo, we obtain in place of the final equation of (2.17), 

(vB)' - (Iu)' = 0, 
or (vB)' — (Iw)' = 0, yielding after integration the relation 

+Iw 



(2.21) 



B 



Substituting in (2.18)-(2.20), we obtain a reduced, planar, ODE in (v,w): 



(2.22) 



(2/i + rj)v' = v{v - 1) + v{p - p_) + ^—((B_ + /w) 2 - v 2 B 2 _), 

2fi v 

fiw' = vw (B- (1 — v) + ito) • 

Mo 
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2.3. The profile ODE as generalized gradient flow. We now recall the general fact [Gj ICS1[ 
CS2] IFR1| [FR2] concerning a hyperbolic-parabolic conservation law 

(2.23) U t + F(U) X = (B(U)U X ) X , U = (%) , F = ( F j) , B - { ° ° 



det£>22 7^ 0, possessing a convex entropy /entropy flux pair 

T) : d 2 r] > 0; q : dq = drjdF 
that is viscosity-compatible in the sense that 

(2.24) ^e{d 2 rjB) > 0, 
that the associated traveling wave ODE 

(2.25) B{U)U' = F(U) - F(U-) - s(U -U-) 
may be written always in the form of a generalized gradient flow 

(2.26) d 2 7]B(U)U' = Vu0{U) 
serving to increase 4>(U) in the direction of positive x. Here, 

(2.27) <t>{U) := sr) - q + dr}{F(U) - F(U-) - s(U -£/_)), 
so that 

(2.28) Vu<t>(U) = d 2 rj(U) (F(U) - F(U-) - s{U - U-)) = d 2 r]B(U)U' 

by direct computation. Likewise, 4>(U)' = Vu4> ■ U' = (d 2 r]B(U))U' ■ U' > by fl2.24| ). We will refer 
to potential <p as the relative entropy (more properly speaking, entropy production). 

Remark 2.1. Evidently, rest points F(U) — F{U-) — s(U — U-) = of the traveling-wave ODE 
correspond to critical points Vi/(p = of the relative entropy. At a rest point the Hessian is 
given by 

(2.29) V^(C/) = d 2 7]{U) (dF{U*) - si) , 

so that, in particular, sgndet = sgn det(dF — si). This gives a connection between the number 
of positive and negative characteristics aj £ o~{dF — si) and the type of the critical point of </>. 



2.3.1. Reduced gradient flow. Using the assumed structure that B has constant left kernel, (2.23), 
as is often the case in applications, we may make a further simplification by the use of entropy 
coordinates. Introducing the entropy variable 

(2.30) W (U) = (Z 1 )(U):=V u n(U), 



globally invertible, by d 2 rj > 0, and noting that dW/dU = d 2 rj, we obtain from (2.26) the more 
useful version 

'0 N 



(2.31) BW = V w <t>, B = B(d 2 rjY 



b 



(2.32) Vw4>={F(U)-F(U-)-a(U-U-)), 

where the block-diagonal form of B follows from vanishing of the first row (inherited from left factor 
B) and compatibility assumption (2.24), which implies also ?R.B > 0. 

Make, finally, the standard assumption (see, e.g., |MaZ3, IZT| IZ2]) that relation 

(2.33) F 1 (U)-F 1 (U-)-s(U 1 -U 1 „) = 

7 



coming from the traveling-wave ODE may be solved for W\ as a function W\ = W(W2) of W 2 , i.e. 
(2.34) det (dw^UiW)) - si) 0. 

Then, defining the reduced potential 



(2.35) 



0(W 2 ) := <ft(W(W 2 ),W 2 ), 



and noting from (2.32) that = for W\ = W(W2), we obtain the relations 



(2.36) 

and 

(2.37) 



V W J = V W2 cP = (F 2 (U) - F 2 {U-) - s{U 2 - U 2 -) 



bW' 2 = V W2 <A, > 



expressing (2.25) as a reduced generalized gradient flow in the parabolic entropy coordinates W 2 
alone. We remark in passing that this implies that (ft = (ft is strictly increasing in positive x 
and not only nondecreasing as shown above. Moreover, we may find (ft directly from (2.36) without 
computing either the full potential (ft or the entropy flux q, which in practice is a great simplification. 

Remark 2.2. In particular, if there is a viscous profile connecting U±, we have 

(2.38) <KU+) = 4>{U+) > $(U-) = 4>{U-). 

2.3.2. Application to MHD. For MHD, we have a viscosity-compatible convex entropy 



/■+00 r+00 

/ p(z)dz+\u\ 2 /2 + v\B\ 2 /2fj, = p{z)dz + \u\ 2 /2 + \vB\ 2 /2^v 

J V J V 



associated with entropy variables ( — p — \ B \ 2 /2fiQ , ui , u 2 , U3 , B\ / hq , B 2 //xo , -B3 j \x§ ) and denoting by 
u = (ui, u 2 ,us) and B = (B\, B 2 , B3) as in ( |2.1[ )p| The associated flux is F(U) = (—u, —Iw,p(v) + 



IB^ 



,,,, . !1q ). The entropy variable is thus 

W := du-n = (-p - \B\ 2 /2/j, , u, w, B/fi ), 



of which the parabolic coordinates are (u,w, B//jq), exactly the ones appearing already in (2.18)- 
(2.20). The corresponding flux density, though we do not need it, is q = \B\ 2 u/ fio + IBw/fiQ + pu. 
Substituting into (2.36) the relation v = u + 1 obtained by integrating the u-equation in the 

^ p{y ) 



traveling-wave ODE, we obtain 
(2.39) V 



V,W,B/fJ,Q ( 



V 



u,w,B / ^io ( 



B 2 -B 2 _ 

V- + 9„„ + v - 1 



2fi 
I(B-B-) 



\ 



MO 

vB-B. 



+ w 
- Iw 



which readily yields 

(ft(u, w, B) 

(2.40) 



p(z)dz — p-(v — 1) + - (w 2 + (v — l) 2 + v 



J 



B 2 -B 2 _ 
Mo 



(B-B)w 



BB_ 



Mo 



Mo 



One checks that d v (ft = d u <ft = p{v) — p- + u + 



w-MB-BJ). 

MO V ' 



B 2 -B 2 _ 
2^o 



, 3 



/' 



(vB — Iw — -B-), and d u 



1 See |Kaw] for related computations in the nonisentropic case. 



2.3.3. The case a = +00. Substituting into (2.36) (2.40) the relation v 
integrating the u-equation in the traveling-wave ODE, and the relation B 

integrating the U-B-equation, we obtain, denoting (f>( 1 •— ! u rl " 

I p(v) — p_ + 



u + 1 obtained by 
- obtained by 



u, w) := <pyu, w 

(B-+Iw) 2 -v 2 B 2 _ 



(2.41) V t 

which readily yields 

(2.42) <f>(v,w) ~- 



I(B-+Iw-B-v) 



+ V - 1 



+ W 



p(z)dz-p^v-l) + l -(v-l) 2 + ^[\ 



flQV 



+ ■ 



Ho 



#2 , l 



Alternatively, this may be obtained from the definition, substituting into (2.40 ) the value B - +Iw for 
B; however, we wish to point out the simplification afforded by working with the reduced problem, 
that is, to emphasize that one need not solve for <p in order to find 0, or for <f> in order to find <j>. 

Remark 2.3. Note that in the above we did not need to compute q or even 77, but only to know the 
entropy variable W, in order to determine the reduced potential by integration of (2.36). Likewise, 

on 

B 2 



computing the full potential <f> by integration of (2.32), we obtain 

B 2 



(v, u, w, B) = u[p{v) — p(l) + 



2fl 



, . (B-BJ) 2 u 2 + w 2 I(B-BJ) 
+ r(v) + ^ — '- + = i J -w, 



2/x 



/'o 



where r satisfies r'(v) = p'(v)(l — v), hence (r(v) — J." p(z)dz)' = r' — p = ((1 — v)p(v))', or 

B 2 -B 2 



(v, u, w, B) = u\p{v) — p(l) + 



2fM 



(B - B-) 2 

+ — + 







u 2 + w 2 



+ 



p(z)dz) + (1 - v)p(v)) 



-w 



I(B - B„ 
2 no 

in agreement up to constant of integration with the formula obtained by direct substitution of r/ 
and q into (2.27). A further substitution yields 

Ml J 

<fi(u + l,u,w, B) = u(p(v) — p(l)) + / p(z)dz — up(v) H — ((v — l) 2 + w 2 ) 



B 2 B 2 
2 Ho 2fi 



BB- 



I 



(B-B-)w, 



Ho Ho 

directly verifying the relation <f>(u + l,u, w, B) = (p(u, w, B). 

2.4. Types of shocks vs. connections. Consider a general system of conservation laws 

U t + F(U) X = {B(U)U X ) X , U G R n 



as in (2.23). Inviscid shock waves correspond to triples (U-, U+, s) satisfying the Rankine-Hugoniot 
conditions 

(2.43) [F(U)] - s[U] = 0, 

where [h] := h(U+) — h(U-) denotes the jump in quantity h across the shock. The type of the shock 
wave is defined by the degree of compressivity 

(2.44) I := dimW(dF([/_) - si) + dimS{dF(U+) - si) - n, 

measuring the number of incoming characteristic modes relative to the shock, where U(M) and 
S(M) denote unstable and stable subspaces of a matrix M, with I = 1 corresponding to the 
classical Lax type, i > 1 nonclassical overcompressive type, and £ < corresponding to nonclassical 
undercompressive type. See |ZHl IMaZ3[ IZ1| for further discussion. 
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At a slightly more detailed level, we define a j-k shock as a shock for which j = n — i_ + 1 
and k = i+ are the indices of the largest positive characteristic speed aj at U— and the smallest 
negative characteristic speed at U+, where a% < ■ ■ ■ < a n denote the eigenvalues of dF(U). Lax 
shocks are associated with a single characteristic family j = k, and we refer to them simply as Lax 
/c-shocks. For overcompressive shocks, j < k, and for undercompressive shocks, j > k, with the 
degree of compressivity i = k — j measuring the difference between j and k. 

Now suppose (as in the present case) that B has a constant left kernel and constant rank, without 
loss of generality 

(2.45) B = (? ? J with b 2 nonsingular, 



y h b 2/ 

and that, if we denote by A the Jacobian matrix of the flux F, 

(2.46) A* := An ~ A^bib^ 1 is nonsingular with real eigenvalues. 

(In the case that there exist a viscosity compatible convex entropy, it may be checked |MaZ3 that 



A* necessarily has real eigenvalues, and det^4* 7^ is equivalent to (2.34).) It follows that traveling 



wave ODE (2.25) can be expressed as a nondegenerate reduced ODE on a manifold of dimension 



r := dimC/2; in the case that there exist a compatible convex entropy, it can simply be expressed 



as the reduced ODE (2.37) in W 2 . 



Suppose further that the shock is noncharacteristic, 

(2.47) det(dF(U± - s)) ^ 0, 
and the endstates satisfy the dissipativity condition 

-OP 2 

(2.48) Ra( - i£A - t 2 B) ± < j^2> > °> for a11 £ ^ R, 

where, here and elsewhere, a{M) denotes spectrum of a matrix or linearized operator M. In the 



case that there exist a viscosity-compatible convex entropy in the vicinity of U±, (2.48 ) is equivalent 
to the genuine coupling condition of Kawashima [Kawj that no eigenvector of A + lie in the kernel 
of B+ , and likewise for A- and £>_ . 

All of these assumptions are satisfied quite generally in applications, in particular for the equa- 
tions of isentropic or nonisentropic MHD with ideal pressure law. See [MaZ4j IZ1| [GMWZ1, 
GMWZ2 for further discussion and examples. 



Lemma 2.4 ( |MaZ3j ). Under the standard assumptions ( |2.45[ ), ( |2.46[ ) ; ( |2.47[ ) ; ( |2.48[ ), U± are 

hyperbolic rest points of the reduced traveling-wave ODE, i.e., have stable and unstable but no center 
manifolds. In particular, for F,B S C , traveling-wave solutions exhibit exponential convergence 

(2.49) \U[x) - U±\ < Ce~ e \ x \ 9 > 0. for x ^ 0. 

Proof. Block matrix reduction and standard invariant manifold theory; see Appendix A, |MaZ3 . 

□ 

Denoting by d+ the dimension of the stable manifold of the r-dimensional reduced ODE at the 
rest point corresponding to U+ and by d- the dimension of the unstable manifold at the rest point 
corresponding to U-, define the connection number 

(2.50) d:=d + +d--r 

measuring the type of the potential connection between rest points U± as a connecting orbit of the 
reduced ODE. Then, we have the following fundamental relation, generalizing the corresponding 
observation of |MP] in the strictly parabolic case. 
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Lemma 2.5 ( |MaZ3] ). Under the standard assumptions ( |2.45[ ), ( |2.46[ ) ; ( |2.47[ ) ; ( |2.48[ ) ; 

(2.51) ^ = d. 
More precisely, 

(2.52) i + = cZ+ + dim 5(A), i_ = d_ + dimW(A)- 



Proof. Results (2.51 )-( 2.52 ) are obtained in |MaZ3j under the additional assumption that there 



exist a connecting profile. However, the proof uses existence only to conclude via homotopy that 
the number of positive eigenvalues dimW(A) of A* is the same at U+ as at [/_, under the weaker 
assumption that A* have real nonvanishing eigenvalues only along the profile. Under our global 
assumption on A*, we have the same conclusions also in the absence of a profile. □ 

That is, the type of the inviscid shock wave determines the type of the potential connection. In 
the simple, planar setting (2.22) of the a = oo case, we have the simple relation that Lax shocks 



correspond to saddle-node connections, overcompressive shocks to repellor-attractor connections, 
and undercompressive shocks to saddle-saddle connections. 

An important consequence is that viscous profiles associated with Lax or undercompressive 
shocks are generically unique up to translation, while profiles associated with overcompressive 
shocks generically appear as part of an ^-parameter family (counting translations). 

2.4.1. Type and orientation. We point out in passing a similar reduction principle at the level of the 
Rankine-Hugoniot equations, this time measuring the parity of d±, or equivalently the orientation 
sgndet(d/([/±) — si) of roots U± of the Rankine-Hugoniot relations. It is sometimes the case that 
certain of the Rankine-Hugoniot equations (2.43) can be solved for certain variables in terms of 
others, that is, without loss of generality, after relabeling U = (U a , Uf,), F = (F a , Ff,), that du b Fb — I 
is invertible, so that Fb(U a ,^(U a )) = 0. In this case, (2.43) reduces to 



(2.53) = F a {U a ) := F a (U a , V(U a )). 

In the present case of isentropic MHD, we will reduce to a scalar equation in the specific volume v. 

Evidently, we have in this case det(dF(U) — I) = det(dbdFb(U) — /) det(dF&(£/) — I), whence, 
since det(dbdFb(U) — I) is real and nonvanishing by assumption, 

(2.54) sgn(det(dF(£/) - /)) = wsgn(det(dF 6 (l7) - I)), u = ±1. 

That is, the orientation of zeros of the full Rankine-Hugoniot relations is determined by the orien- 
tation of zeros of the reduced Rankine-Hugoniot relations (2.53). We make use of this later to help 



determine the types of rest points by consideration of a scalar reduced relation. Similar reasoning 
is used in |FR1] for a planar reduced relation, looking at orientations of intersections of nullcline 
curves (equivalent to orientation of zeros of the planar reduced condition). 



See Remark 2.1 and Appendix |A| for related observations. 



2.5. The Evans function and stability. We conclude these preliminaries by a brief discussion 
of stability of general traveling-wave profiles, as determined by an Evans function, or "generalized 
spectral stability" condition. Throughout this section, we make the general assumptions (2.45), 



( [2746] ), fl2747| > ] fl2748| ) of |MaZ3l IMaZ4l [ZT] . as hold in particular for the MHD equations studied 
here. We add to these the further assumption of symmetric- dissipative hyperbolic-parabolic form 



(S) There exist coordinates W for which (2.23) becomes G(W) t + F(W) X = (B(W)W, 



with dG symmetric positive definite and block-diagonal, dF\\ symmetric and either negative or 
positive definite, and B block-diagonal, with 5R^22 positive definite. (Here and elsewhere, $l(M) := 
(1/2) (M + M*) denotes the symmetric part of a matrix or linear operator M.) 
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This structure guarantees the minimal properties needed to carry out an analysis, in particular 
that the nonlinear equations be local well-posed and that the linearized equations generate a C° 
semigroup; see |Z2| GMWZ1, GMWZ2 for further discussion. It is implied by existence of a 
viscosity-compatible convex entropy together with the condition that a (A*), real and nonzero by 
assumption, be also strictly positive or strictly negative, a minimal further requirement since A* 
in most applications is a scalar multiple of the identity. In particular, (S) and all other hypotheses 
are satisfied for the equations of MHD with ideal gas equation of state [MaZ4, Zl under the single 



condition (2.47) of noncharacteristicity. 

Linearizing about a stationary wave U = U(x) of ( 2.23| ) (stationarity may always be achieved 



by a change to coordinates moving with the wave), we obtain linearized evolution equations 
(2.55) U t = CU := (BU X ) X - (AU) X , 

where A and B depend on x, converging asymptotically to values A(±oo) = dF(U±), £>(±oo) 



B(U±). By asymptotic convergence (2.49) and dissipativity, (2.48), we find from a standard result 
of Henry |He] equating essential spectrum of asymptotically constant coefficient operators to that 
of their limiting constant-coefficient operators, that 

a ess {C) C {A : 5RA < -#|3A| 2 } 

for some 9 > 0, where a ess denotes essential spectrum (defined as the part of the spectrum not 
consisting of eigenvalues); see |AGJ[ IGZ[ IZlj . Moreover, this bound is sharp; in particular, A = 
is in the limit of the essential spectrum. At the same time, A = is always an eigenvalue of C, due 



to translational invariance of the underlying equations (2.23), with associated eigenfunction U'. 

The fact that there is no gap between the spectrum of C and the imaginary axis makes this a 
degenerate case for which linearized and nonlinear stability analysis is trickier than usual. In the 
standard case of a sectorial operator for which there exists a spectral gap, one may conclude bounded 
linear stability from the spectral stability conditions of (i) nonexistence of unstable eigenvalues 
5RA > 0, and (ii) semisimplicity of neutral eigenvalues 3ftA = 0; indeed, these are necessary and 
sufficient. However, here we have a nonsectorial operator with no spectral gap. Moreover, for the 
eigenvalue A = embedded in the essential spectrum of C, it is not clear even what is the meaning 
of semisimplicity; see discussions of |ZH| IMaZ3| IZ1| IZ2j . 

Nonetheless, as shown in [GZ] IZH| IMaZ31 IMaZ4| IZlj . one can extract a simple necessary and 
sufficient condition for stability analogous to (i)-(ii) in terms of the Evans function -D(A) associated 
with C, a Wronskian 

(2.56) D(X) := det(Wf, . . . , W k , W+ +1 , W N )\ X=0 

defined in terms of analytically-chosen bases {W-[~, . . . , Wj^}(\, x) and {VF fc t +1 , . . . , Wn)}(\, x) of 
the manifolds of solutions decaying as x — > oo and x — * +oo of the eigenvalue equations (£— X)w = 
written as a first-order system 

(2.57) W' = A(x,X)W, 

where W is an augmented "phase variable" including w and suitable derivatives. By standard 
considerations, this may be defined on the complement of a ess (C); a more detailed look shows that 
D permits an analytic extension to the boundary of this set- in particular, to the nonstable complex 
half-plane {3?A > 0}. For details of this construction, see, e.g., |AG J| iGZl IZ1 | lHuZ2j : we give some 
further discussion also in Section [5] and Appendix [D] 

Evidently, away from the essential spectrum o~ ess (£), the Evans function vanishes at A if and 
only if A is an eigenvalue of C, corresponding to existence of a solution of the eigenvalue equations 
decaying at both x — > ±oo. Indeed, the multiplicity of the root is equal to the multiplicity of the 
eigenvalue |GJ1[ IGJ2[ IMaZ3[ IZ1] . The meaning of the multiplicity of the root of D at embedded 
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eigenvalue A = is less obvious, but is always greater than or equal to the order of the embedded 
eigenvalue |MaZ3llZl] . 

By the discussion in Section 



2.4 



in particular relation (2.51), a traveling-wave profile U lies in 



an ^-parameter family of nearby solutions, where £ is (by dimensionality) at least min{l,£}, where 



I is the degree of compressivity defined in (2.44), with equality in the case that the connection is a 



maximally transversal intersection of the unstable manifold at U- with the stable manifold at U+. 
Assume for simplicity the typical case that equality holds, 

(2.58) £ = mm{l,£}, 

and the manifold of nearby solutions is smooth. Then, the stability condition is 

(D) D has precisely £ roots on the nonstable half-plane {A : KA > 0}, necessarily at A = 0. 

This is analogous to the stability condition in the standard sectorial case, with nonvanishing away 
from A = corresponding to the standard spectral condition (i), and vanishing to order £ at A = 
indicating that the multiplicity of this zero is accounted for entirely by genuine eigenfunctions 
corresponding to variations of the traveling wave connection along the ^-parameter family of nearby 
solutions, a generalized version of semi-simplicity |ZH, MaZ3l [ZT] . 

2.5.1. Linear and nonlinear stability. We have the following basic results relating the Evans con- 
dition (D) to stability. 



Proposition 2.6 (|MaZ3j). Under the standard assumptions (2.45), (2.46), (2.47), (2.48), and 
assuming an £-parameter family of traveling -wave solutions near U, (D) is necessary and sufficient 
for linearized stability from L 1 D L°° — ► LP of U, all 1 < p < oo, defined as \e Lt f\hv < Cl/liini"- 



Proposition 2.7 (|MaZ4, RZ ). Under assumptions (2.45)^(2.48), (S), and definining £ as in 
(2.58), the Evans condition (D) implies, first, existence of a C 1 family of nearby solutions {u a }, 
a G R , and, second, nonlinear time- asymptotic orbital stability, in the following sense: For any 
solution U of (2.23) with initial difference Eq := ||(1 + |x| 2 ) 3//4 (f/(-, 0) — f7) || j^s sufficiently small 



and some uniform C > 0, U exists for all t > 0, with 

(2.59) ||(1 + \x\ 2 f /4 {U{-,t) - U{- - st))\\ H B < CE (stability). 
Moreover, there exist a{t), such that 

(2.60) \\U(; t) - U a( t\- - st))\\ LP < CE (1 + t )-(i/2)(i-i/ P ) ) 
and 

(2.61) \a(t) -Oool < CEoil + ty 1 / 2 , \a(t)\ < CE (1 + 
for all 1 < p < oo. (phase-asymptotic orbital stability). 

A similar result holds in the mixed, under-overcompressive case that the family of nearby trav- 
eling waves has dimension different from (necessarily greater than) £; see [RZ . 



2.5.2. The integrated Evans condition. Noting that C = d x {Bdx — A) is in divergence form, we 
may conclude for any A 7^ that satisfaction of the eigenvalue ODE (£ — X)w = for an solution w 
decaying exponentially in x up to one derivative implies that w{x) := w(y)dy is also bounded 
and exponentially decaying, and satisfies the integrated eigenvalue equation 

(2.62) (C-X)w = 0, 

where C := Bd 2 — Ad x . Associated with C is an integrated Evans function D(X), which like D 
may be defined analytically on the nonstable half-plane {A : 5RA > 0}. This permits the following 
simplified stability condition, in practice easier to verify. 
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Proposition 2.8 (|ZH, MaZ3j). Under assumptions (2.45)^(2.48), in the Lax or overcompressive 
case, the Evans condition (D) is equivalent to the integrated Evans condition 



(D) D is nonvanishing on the nonstable half-plane {A : 3?A > 0}, 
and in the undercompressive case to 

(D') D has on the nonstable half-plane {A : 5RA > 0} a single zero of multiplicity one at A = 0. 

In the Lax and overcompressive cases that are the main focus of our investigation here, the 
change to integrated coordinates has the effect of removing the zeros of D at the origin, making 
the Evans function easier to compute numerically and the Evans condition easier to verify. 



3. Rankine-Hugoniot Conditions 

The Rankine-Hugoniot conditions for isentropic MHD are, in the notation u = u\, B 

w = (u 2 ,u 3 ), 



(S 2 ,S 3 ) 



(3.1) 
(3.2) 

(3.3) 
(3.4) 



-s\u 



-s\w\ 



= w 






p + 




' B~ 


= I 






/-to. 


= I\w}. 



B*_ 
2/xo 



Under the scaling (|2.15|), we have s 

-- 0, w 



invariance), we may take u 
rotational invariance) ■ 
which, so long as 



1 , t>_ = 1 , and without loss of generality (by translation 
= 0. Last, we may take without loss of generality (by 
0, whereupon we obtain from (3.3)-(3.4) that [S3] = [vB 3 ] 



0, 



(3.5) v + + = 1, 
gives, finally, 

(3.6) S 3 . 
Collecting, we have the normalizations 

(3.7) s = —1, i>_ = 1, it_ = 0, W2 



B 



3+ 



0. 



0, w 3 _ = w 3+ = 0, S 3 



B 



3+ 



0. 



To generate all possible shock profiles, up to invariances of the equations, we shall vary /, S2+, 
without loss of generality nonnegative, and v+, without loss of generality between and 1 (since we 
can always arrange that v— = 1 correspond to the rest point with larger v value), and solve for the 
remaining coordinates u+, B2-, and the parameter a appearing in the pressure law. Parameters 
that will be important in the whole study are 



(3.8) 



J 



(S 2 



2/u 



and K :- 



J 



Mo' 



Bl 



2v- 



'M0 



{I? 



{I? 



(Note that, under the rescaling that we used, I = 
in the original coordinates.) 

Remark 3.1. In the excluded case u+ = V- = 1, profiles are prohibited by entropy consideration, 
(2381). 
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Remark 3.2. Note that it does not follow in general that B3 = or 103 = 0, but does follow 
when profiles are unique, i.e., in the Lax or undercompressive case, and one such profile is known 
to exist. (Recall the discussion of types of shocks and relation to uniqueness of profiles in Section 



2.4) 



Proposition 3.3. Under the normalizations (3.7), for each < w+ < 1 and I, B 2+ > 0, the 



Rankine-Hugoniot equations (3.1) -(3.4) have a unique solution 

K\„ K 

T 



(3.9) 



Hi 



1, B 2 



1-K 



B 2 



l-v+ 
1 - K 



£2, 



l-v+ 



B\ + (1 + v+ - 2K) \ / 1 - v + 



This is physically meaningful if and only if a > 0, or 
(3.11) 



1 - J 



(l + v+-2K) \ 
{v+-K) 2 )' 



l<v + -K2{K-l) + {l-Kf(^f). 

2+ 



(For K > 1/2, this gives no restriction. For K < 1/2, B\ < 



2 M0 (1-E-) 2 
\-2K 



or J < 



1-2K v 



Proof. From [u] = [v], we obtain immediately u + = v+ — 1. Expanding [Bv] = I[w] = K[B2] and 
solving, we obtain 

-v + -K- 



Bo 



1 - K 



Bo 



From [w] = (K/I) [B2] , we then obtain w + 



T ( i-x ) • Finally, from the remaining condition 



p + 



2fi 



, we obtain 



(3.12) 



\ p ] = (l-v + )-(l/2n )(Bi + -B 2 2 _), 



yielding (3.10) and (3.11). 



□ 



Remark 3.4. So far, we have made no restriction on dimension or a, so our analysis of the 
Rankine-Hugoniot conditions holds for the general three-dimensional isentropic case. 



3.1. Global rest point configuration. Proposition 3.3 gives a convenient means for stepping 



through the possible shock connections, and is the main method we will use to generate shocks in 
our numerical investigations of shock stability. For the study of the existence problem it is more 
useful to take a global point of view, fixing a left state and speed in the unrescaled coordinates, 
and studying the configuration of rest points (possible right states) in the resulting traveling-wave 
ODE. In the rescaled coordinates, this amounts to fixing I, B 2 ~, and a, or, equivalently, the more 
convenient parameters (J,K,a), and solving for all possible v+. 

Proposition 3.5. In the parallel case J = 0, for < a / 7 -1 and < K 7^ 1, there exists a unique 
parallel solution 7^ 1 satisfying = g{v) := p{v) — p(l) + v — 1, with associated magnetic field 
B 2 * = 0. If K is not between v* and 1, then these are the only rest points, with t>* corresponding to 
a saddle and 1 to a repellor if K < < 1 and 1 corresponding to a saddle and v* to an attractor 
if < 1 < K. If K lies between v* and 1, then 1 corresponds to a repellor and v* to an attractor 
and there are two additional nonparallel saddle-type rest points 



v = K, u = K-l, B = ±^2fx {p(K) - p(l) + K-1), w = KB/I. 
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Proof. We have g(l) = and g'(l) = p'(v) + 1 = —7a + 1 7^ by assumption. Since g is evidently 
convex, and g — ► +00 as v — ► 0, +00, we find that there is precisely one other root 7^ 1. There are 
a further two solutions v = K w = —KB2/I, B 2 /2^lq = — ([p] + [f]) = —g(K), which are physically 
relevant only if g{K) < 0, or (by convexity) K lies between 1 and v*. The types of the rest points 
may be obtained by straightforward computation |BHZ]. □ 

Proposition 3.6. In the nonparallel case J > 0, for a > and < K 7^ 1, rest points of 
traveling-wave ODE (2.13), or, equivalently, right states satisfying the Rankine-Hugoniot equations 
(3.1)-(3.4) with s = —1 and V- = 1, correspond to roots v = v+ of 



(3.13) 



f(v) := p(v) - p(l) + J 



(v - Kf 



of which there are at most two greater than K and at most two less than K. For all except a measure- 
zero set of parameters, there are exactly two or four roots in total, consisting of an attractor v\ and 
a saddle V2 ordered as v 1 < V2 < K, a saddle V3 and a repellor v<± ordered as K < v% < V4, or both, 
with (u,w,B) values determined by (3.9). Moreover, the relative entropy 4>(vj) decreases with j. 



Proof. Combining ( |3.9| )(ii) and (3.12), we obtain (3.13). Noting that /, since p is convex, is convex 
on (0,-ff) and (K, +00), with f{v) — > +00 as v — ► 0, K, +00, we find that / can have at most 
two roots on each of the intervals (0,K) and (K, +00). Noting that / is monotone increasing in 
a on (0, K) we find for each fixed (J, K) that there are at most two values of a for which / has a 
double root, hence, for all except this measure zero set of parameters, there are exactly two or four. 



Applying the reduced orientation principle (2.54) together with the reduced type relation (2.52) 



we find using the fact that /' changes sign between two roots on one side of K that one must be of 
saddle type and the other of node type. 

Finally, tracking down the orientations of intermediate transformations, which change sign as v 

crosses K, by the relation w = f^p^)^, or, more simply, directly computing the sign of the 



determinant of the 2x2 matrix arising from the linearization of the planar ODE (2.22 ) about the 
rest points (v,w), we find that the largest root > K and the smallest root < K are nodes, and the 
others saddles. Computing the trace of the 2x2 coefficient matrix of the linearized system, we 
find that the largest root > K is a repellor and the smallest root < K an attractor. 
Alternatively, and much more simply, recalling the formula (2.41) for V vw (j), solving 



(3.14) 



= V, 



/(£>_ + Iw - B_v) 



+ w 



w(v - K) LB_(1 



for w — 



„ u,-k~\ ; an d substituting into V v <p, we find after a brief computation that, along this 

nullcline, d<p/dv = f(v), hence the relative entropy (j) is decreasing with respect to v between rest 
points lying on the same side of K, again identifying nodes > K as repellors and nodes < K as 
attractors for the flow of the planar traveling-wave ODE. (Recall that increases along the flow, 
with rest points of the flow corresponding to critical point of <fi.) Finally, taking without loss of 
generality B- < 0, note that, by (3.14), at v = K, V w <fi = — IB ~^~'"^ > for all w, so that the 

limiting value of q> as v — * K + on the negative-u> nullcline branch for v > K is less than the limiting 
value of as v — * K~ on the positive-w; nullcline branch for v < K, verifying decrease with j of 



for all j and completing the proof. 



□ 



Remark 3.7. Note that the above argument depends only on the general properties of the pressure 
law p of convexity, blowup at at rate at least 1/v and decay as v — ► +00, and not on the specific 
form of a polytropic gas law, hence our conclusions extend to general pressure laws of this type. 
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Remark 3.8. The parallel and nonparallel cases can be combined, associating rest points to roots 
of the continuous function f(v) := (v — K) 2 f(v). In all cases, there is a two-parameter bifurcation 
at K = 1, with three rest points collapsing at v = K = 1. 

Remark 3.9. Though we carried out our analyis for the planar system arising through the choice 
a = oo and the restriction to two dimensions, our conclusions on the number and type of states 
satisfying the Rankine-Hugoniot conditions apply to the general three-dimensional isentropic case. 



That is, through the relation (2.52) we are able to make quite general conclusions on types of 
shocks by examination of the simple planar realization of the two-dimensional traveling-wave ODE. 
Indeed, our final argument determining the type of rest points for the planar system by looking 
along nullclines of amounts to a further reduction to the scalar realization obtained by setting 
fx = as well as a = oo. 



Factoring out the root v = 1, we may examine instead roots of 

/>! = Pfo)-P(l) J(l+v-2 K) 
1 v-l {v-Kf 



(3-15) f(v) := fH = —. ?' + "J ^2 ' + !" 



Remark 3.10. In describing the possible four rest point configurations in the nonparallel case 
J > 0, we may (by rescaling if necessary) without loss of generality consider only the case that 
v_ = 1 is the largest rest point: that is, /'(l) > and there is a rest point v+ < K < 1. Fixing u + 



and K, and letting J vary, we obtain by (3.10) that 



, 1 - v+ \ , / 1 - v+ \ (1 + v+ - 2K) 

(3.16) < a = c - dJ, where c = — + > and d - ' 1 



\) V^-T _ i7 (v + - K) 

Thus, since K < 1, 

2 J / 2 

(3.17) < /'(l) = - 7 a + 1 - — - = (1 - 7 c) + J(7d - 



implies either d > 0, in which case J < c/d = ^^^^ by (3.16), or else d < 0, in which case 

J < _ 7t j_| 1 2/(i-_g) • ^ e same considerations hold whenever there exist Lax 1-shocks, or, without 
loss of generality (by rescaling the largest rest point to value = 1) K < 1. That is, it is sufficient 
to consider a bounded parameter range (a, J) in studying four rest point configurations or Lax 1- 
shocks, for K bounded away from 1. This is important for numerical explorations, in which the 
parameter range is necessarily finite. 

3.2. Four rest-point configurations. To aid our later numerical investigations, we give a simple 
description of the set of parameters (J, K, a) for which four-rest point configurations appear, and 
with them the possibility of intermediate, over compressive, and undercompressive shocks, without 
loss of generality taking K < 1 by rescaling if necessary so that V- = 1 is the largest root of /. 

Proposition 3.11. For < K < 1, the set of J > and a > for which there exist four solutions 



of the Rankine-Hugoniot equations (3.1 )-( |3.4[ ) (equivalently, four rest points of traveling-wave ODE 
( |2.18 )-(2.20) or (2.22)), except for the measure-zero set of values a = a*(J,K) := for which 

/'(l) = ; consists of a connected set 

(3.18) {(J, K) G n := Tlx U K 2 , < a < A(J, K)} 
for some A(J, K) > 0, where 

(3.19) Tlx ■= |0 < K < -, < J < 2R _ - } and K 2 := < K < 1, < jj. 

17 



Proof. As there is always a rest point with V- = 1, existence of four rest points is equivalent (except 
on the measure zero set of parameters for which degenerate roots appear), to existence of a second 
rest point with v £ (0,K), i.e., a root of /. Since / is monotone increasing in a for < v < K, 
this consists of an open interval a £ (0, A(J, K)), with A > only if / has a root in (0, K) for the 
limiting values (J,K,0). 

Multiplying / by (v — K) 2 /(l — v) reduces this question to existence of a root v £ (0,K) of 
the quadratic q(v) = — (v — K) 2 + J(l + v — 2K), q'(v) = —2{v — K) + J. We readily compute 
that q(K) = J(l - if) > for J > and g'(if) = J > and that q(6) = -K 2 + J(l - 2K) and 
g'(0) = 2K + J > 0, so that the only way there can be a root of q on < v < K is if q(0) < 0, 
or J(l — 2K) < K 2 . The set of (J,K) satisfying this condition is easily seen to correspond to the 
connected set IZi U 7^2 • D 

Remark 3.12. In the case v > K > 1, / is monotone decreasing with a, and so we find that the set 
of parameters generating four rest point configurations (ignoring the measure-zero set corresponding 
to /'(l) = 0) is, rather, of form a > A(J,K) for arbitrary J > 0, K > 1. The set of (J, if) for 
which four point configurations appear for all a > 0, i.e., A = 0, is readily seen to be 1Z% := 
{K > 1 and J > 4(if - 1)}. For, in this case, q{K) = J(l - K) < but q'{K) = J > 0. 
Meanwhile, g(+oo) < as well, so that the only chance for a root v > K is that the maximum 
value of q be positive. Solving q'{v) = for the critical point v m ax = J/2 + K, we find that 
q(v ma x) = J(J/4 + 1 — K)-, which is positive precisely for J > 4(K — 1). 

Recall, for four rest point configurations, we may take without loss of generality K < 1. 

3.3. Two-dimensional shock types. Restricting to two dimensions, we find that shocks connect- 
ing rest points in decreasing order of v are of Lax 2- type for both values < K, of Lax 1-type for both 
values > K. Shocks connecting the largest v value to the smallest are overcompressive 1-2 type, 
while shocks connecting the largest v-v&lue > K to the largest v-value < K are Lax 1-type and 
shocks connecting the smallest u-value > K to the smallest < K are Lax 2-type. Shocks connecting 
the two middle (nonextremal) values of k are undercompressive 2-1 type. In the terminology of the 
literature (Gj ICS1| ICS21 IFS] , all shocks bridging across the value K are called intermediate shocks] 
as shown above, these may in principle be of Lax, overcompressive, or undercompressive type. 

In our main parameter range (monatomic or diatomic gas with standard viscosity ratio for non- 
magnetic gas), only Lax and overcompressive type appear to have profiles for the two-dimensional 
a = oo case considered here. 



3.4. The three-dimensional case. The results of Propositions 3.3 3.5 and 3.6 extend by rota- 
tion to the full, three-dimensional case to yield the same basic 2-4 rest point configuration, with 
all rest points confined to a rotation of the planar case. The single exception is in the parallel case 
J = 0, for which the data corresponding to u_ = 1, and to the attractor v* < K if it occurs, is 
rotation invariant; in this case, the intermediate rest points extend by rotation to yield a circle of 
rest points. Associated intermediate shocks are called degenerate type |FSJ; in case J / 0, they 
are called nondegenerate type. The double-cone configuration arising from rotation of the four rest 
point parallel configuration, and associated interesting bifurcations, are discussed in |FS]. 

Considered as waves of the full three-dimensional system, shocks connecting rest points in de- 
creasing order of v are of Lax 3-type for both values < K, of Lax 1-type for both values > K. 
Shocks connecting the largest v value to the smallest are overcompressive 1-3 type, while shocks 
connecting the largest v-v&lue > K to the largest v-value < K are overcompressive 1-2 type and 
shocks connecting the smallest v-v&lue > K to the smallest < K are overcompressive 2-3 type. 
Shocks connecting the two middle (nonextremal) values of k, undercompressive when considered 
as two-dimensional waves, are in three dimensions of Lax 2-type, or Alfven waves. Thus, in three 
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Figure 1. Region of four rest point configurations, a = 0, in (K, J) space. 



dimensions, only Lax or overcompressive shocks appear. That is, undercompressivity is an artifact 
of the restriction to two dimensions. 

In our main parameter range (monatomic or diatomic gas with standard viscosity ratio for 
nonmagnetic gas), Lax 2-shocks do not appear to have profiles when resricted to two dimensions, 
for the o~ = oo case considered here. In the full, three-dimensional a = oo case, therefore, Lax 
2-shock profiles if they exist must be nonplanar in the sense that they leave the plane of the rest 
point configuration. Likewise, overcompressive shocks, besides the planar connections studied here, 
admit also nonplanar connections when considered in the full, three-dimensional setting. We shall 
not study such genuinely three-dimensional profiles here, restricting attention to planar profiles that 
can be studied within the two-dimensional framework. For discussion of fully three-dimensional 
phenomena in the related nonisentropic case, see (Gj ICS1[ ICS21 IFS] . 



4. Existence of profiles 

In this section, we describe the possible viscous shock profile connections for the various rest 
point configurations described in Section K3j Typical phase portraits (determined numerically) for 



the two variable system (2.22) with a = oo are graphed in Figures PZl pi andp? 



4.1. The parallel case, J = 0. 

Proposition 4.1. In the parallel case J = 0, for < a / 1 and < K ^ 1, and a = oo, assuming 
without loss of generality that 1 is the largest rest point of the traveling-wave equation, there is 
always a profile connecting V- = 1 and the unique parallel rest point / 1, w* = 0, B^* = 0, 
which is of Lax 1-type if K < v*, Lax 3-type if K > 1, and overcompressive type if v* < K < 1 . 
In the latter case, there are Lax connections from repellor i;_ = 1 to the additional saddle-type rest 
points 

v = K, u = K-l, B = ±^2n (p(K) - p(l) + K-1), w = KB/I, 

and from these saddle-type rest points to the attractor v*, whose orbits bound a four-sided region 
foliated by overcompressive connections. 



Proof. In the parallel case, (2.22) reduces to 

/ J™ 2 

(41) {2^ + r,)v =vh(v) + —, 

fiw' = (v — K)w, 
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Figure 2. Nullclines and phase portrait for typical parallel case, parameters v+ = 
0.1, I = 0.7, B + = 0, and p = r = 1. 



where h{v) := (v — 1) + (p — p_) is convex and vanishing at u#, 1, hence negative for v £ (v*, 1). 
Setting uu = 0, we find that there is a monotone decreasing solution connecting v_ = 1 to u*, which 



has the type described by the results of Proposition 3.5 

As h < 0, the nullclines = ±^ —v 2 h(v) / J for u' are well-defined for f G (u*, 1), bounding a 
lens-shaped set 7£ between u # and 1, passing through the saddle-type rest points at v = K and 
pinching to a single point at v = 1, u#, within which i/ < 0. Noting that sgnu/ = ±sgnu; for v ^ K, 
we find that this region is invariant in backward (resp. forward) x, whence, starting at the saddles 
and integrating in backward (resp. forward) x along the stable (resp. unstable) manifold, we find 
that the orbit remains for all x in 1Z, hence must connect to 1 (resp. v*), verifying existence of the 
bounding Lax-type connections. Starting at any point {K, w) lying on the open interval between 
the two saddles and integrating in both forward and backward x, we likewise find that the orbits 
are all trapped in 1Z for all x, so generate a one-parameter family of overcompressive connections 
filling up TZ. See Figure [2} □ 

4.2. Existence of Lax-type profiles, J > 0. 

Proposition 4.2. In the nonparallel case J > 0, for a > and < K ^ 1, with a = oo, rest 



points Vi < Vj of traveling-wave ODE (2.22) lying on the same side of K always admit a Lax-type 
profile, which, moreover, is monotone in both v and w. 

Proof. Without loss of generality, let the rest points be v+ < 1 and V- = 1. Rewriting (2.22) as 

(2/i + rj)v' = vh(v) + (^-+ Iw )\ 
(4-2) 2/W \ , 

fiw ={v — K)w , 

Mo 

where h(v) := p(v) —p~ + v — 1 — J is convex and (since vh(v) + = 0) negative at v = v+, 1, 

hence negative on (v+, 1), we find that the nullclines Iw = —B_ ± \J — 2fiQV 2 h(v) for v' for v' are 
well-defined for v G (v+, 1). Likewise, the nullclines w = ^t^~w for it;' are well-defined for v on 
either side of K, forming two disconnected branches asymptotic to the line v = K. 

Case K < v+. In this case we find that the rest points v+, 1 must lie on the intersection of the 
lower branch of the nullcline v' = 0, and the righthand (> K) branch of the nullcline w' = 0, and 
these nullcline branches have no other intersection (else there would be a third rest point for v > K, 
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impossible by Proposition 3.6). Looking at asymptotics, we find that the nullcline w' = must 
lie above the nullcline v' = for v G t(v+, 1), with the two curves forming a lens-shaped region 1Z 
between v+ and 1, within which v' < and w' < 0. Looking along the boundaries, we find that 
the vector field (v',w') points out of 1Z, so that 1Z is invariant in backwards x. Thus, integrating 
backward in x from v = v+ along the stable manifold, we find that there exists a connection to 
v = V- = 1, which is monotone decreasing in v and w. 

Case K > 1. A symmetric argument yields existence in case K > 1, again with v monotone 
decreasing and w monotone increasing, this time via invariance in forward x. See Figure [3j □ 

Remark 4.3. The argument above may be recognized as the same one used to prove existence of 
nonisentropic gas-dynamical profiles in |Gi] . It should be possible to obtain this result alternatively 
by a relative entropy argument as in [G] for the nonisentropic case, showing in case K < v+ that 
the level set of </> through v+ encloses V— = 1, yielding existence by a Lyapunov- function argument 
in backward x; this would apply also for a finite. 

4.3. Existence of intermediate shock profiles, J > 0. 

Proposition 4.4. Set r := /i/V. In the nonparallel case J > 0, with a = oo, for each fixed (a, J, K) 
with a > 0, J > 0, and < K ^ 1 for which there exist four rest points v\ < V2 < K < v% < of 



traveling-wave ODE (2.22), there exists a value r* = r*( J, K, a) > such that: (i) for r < r* ; there 
exist no intermediate shock profiles (i.e., the only connections are regular Lax profiles between V2 and 
vi and V4 and V3 as described in Proposition \4 ■ Sty ; (ii) for r = r* ; there exists an undercompressive 
profile connecting V3 to V2, monotone decreasing in v and increasing in w, and no other intermediate 
shock profiles; (Hi) for r > r* , there exist intermediate Lax connections from V3 to v\ and v 4 to v 2 , in 
general not monotone in v or else not montone in w, and a one-parameter family of overcompressive 
profiles from to v\, in general not monotone in v or w, with no other intermediate shock profiles. 



Proof. Referring to Figure |3[b), rewrite (2.22) again as 

, - (B.+Iw) 2 
tv =vh(v)-\ , 

(4-3) 2W \ , 

fJLW = (V — K)W . 

MO 
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r = 2/i + 77, where (see proof of Proposition 4.2 ) := — p_ + v — 1 — J is convex, negative 
on (ui, W4), and goes to +00 as v — > 0, +00. Denote by v < v 1 and v > V4 the two points at which 
h vanishes. 
The nullclines 

Iw = -B- ± ^-2^ v 2 h{y) 

for v' evidently are well-defined on v £ (v,v), together forming a simple closed curve C enclosing 
a region on which v' < 0, as seen in Figure 3(b). Likewise, the nullclines w = ^~t^~^t for w' are 
well-defined for v on either side of K, forming two disconnected branches B\ and B2 asymptotic to 
the line v = K, Figure [3^b). 

The arc formed by the portion of C from the rest point at v% to v together with the portion of 
the v axis from v to 0, the portion of the w axis from to the intersection of the w axis with B\, 
and the portion of B\ between the intersection of the w axis with B\ and the rest point at v\, the 
Lax connection between the rest points at v \ and V2, and the portion of B\ between the rest point 
at t>2 out to w — > +00 form a barrier to the flow in forward x, through which an orbit initiating 
inside C cannot cross, as, likewise, does the arc formed by the Lax shock between the rest points 
at v\ and v<i together with the portion of B2 extending from the rest point at V2 to w — > —00. 

Thus, the orbit O initiating along the unstable manifold of the rest point at V3 pointing in 
decreasing v-w directions, and thus initially lying inside C, must either (a) strike the arc between 
(v, 0), after which, being trapped between B\ and C, it must asymptotically approach the rest point 
at v\] (b) strike the arc of B\ between the rest points at v\ and V2, after which, being trapped 
between this arc and the portion of C between the rest points at vi and V2, it must asymptotically 
approach the rest point at v\\ (c) remain within the interior of C and to the right of B\ for all time, 
asymptotically approaching the rest point at V2\ or, (d) exit the interior of C along the arc between 
the rest point at v% and the rest point at vi, after which it remains trapped outside of C with v 
increasing monotonically to +00. 

Depending whether the orbit approaches the rest point at vi, approaches the rest point at 
V2, or takes v — > +00, we are in cases (iii), (ii), or (i) of the proposition. But, these cases are 
distinguished by the location along the arc C formed by the portion of the upper branch of C 
lying below B\ together with the portion of B\ lying below C at which the orbit O exits the part 
of the interior of C lying below B%, with the locations corresponding to different cases ordered in 
clockwise fashion along C Noting that the signs of v' and w' are constant while O; remains inside 
C (recall that it is trapped to the right of Si), and are given respectively by r _1 and fi~ 1 times the 
righthand sides in ( |4.3[ ) , we find that the exit point moves strictly clockwise along C monotonically 
as r = r -1 //^ 1 = /jl/t increases. 

Thus, as asserted, there is a unique value r = r* for which O exits at the rest point at V2, 
corresponding to an undercompressive connection. For r < r*, O exits to the right of the rest point 
at V2 , going off to infinity, and for r > r* , O exits to the left of the rest point at 1)2 , asymptotically 
approaching the rest point at v\, corresponding to an intermediate Lax connection and case (iii). 

Note, in case (iii), that the existence of this intermediate Lax connection means that, applying 
a symmetric argument in backward x to the orbit originating from the saddle rest point at V2, we 
find that it remains trapped within C for all negative x, approaching asymptotically as x — » —00 
the rest point at v\ The four Lax connections enclose an invariant region, within which all orbits 
must be overcompressive profiles connecting the rest points at v\ and U4. It is clear that in general 
the intermediate Lax profile may leave either the interior of C or the region below B\, hence may 
be nonmonotone in v or w but not both. Similarly, we find that the members of the family of 
overcompressive profiles are in general nonmonotone in v or w (and sometimes both). 
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In case (ii), or case (c) above, the profile remains for all x in a region for which v' < and w' > 0, 
hence the profile is monotone decreasing in v and increasing in w. This completes the description 
of the phase portrait in cases (iii) and (ii), finishing the proof. □ 

Remark 4.5. Evidently, there is nonuniformity in the behavior of r*, in view of fact that the 



parallel case J = is always in case (iii) of Proposition 4.4 by the result of Proposition 3.5 That 
is, r* — ► oo as J — > for fixed K and a. 



4.3.1. Singular perturbation analysis. The results of Proposition 4.4 may be illuminated somewhat 
by formal singular perturbation analyses as r — * and r — ► +oo: equivalently, taking — > with 
r = 1 fixed, or r — > with jx = 1 fixed. In the limit as /i — ► 0, the phase portrait for J > reduces 
to slow flow along the w' = nullcline ,8 (notation of the proof above), with fast flow involving 
jumps in the vertical w direction. We find that regular Lax connections are accomplished by slow 
flow along C, but there are no further intermediate shock connections since the branches B\ and 
£>2 of B are separated by the vertical line v = K. See Figure |4jh). In the special case J = 0, 
the hyperbolae Bj degenerate to the connected union of v = K and w = 0, allowing intermediate 
connections both from the rest point at v± to the rest points at V2 = ^3 and from the rest point at 
Vi to the rest point at v\. 

In the limit as r —* 0, the phase portrait reduces to slow flow along the v' = nullcline C 
(notation of the proof above), with fast flow involving horizontal jumps in v. We find that Lax and 
intermediate Lax connections may all be accomplished by slow flow along C, with fast flow filling 
in the overcompressive family. See Figure |4^a). 

Finally, note that as r goes from + (/z —* limit) to +oo (r —* limit) the relative orientation 
as measured by a Melnikov separation function along an appropriate transversal of the unstable 
manifold pointing to the left at r = + of the rest point associated with t>3 and the stable manifold 
entering from the right at r = + of the rest point associated with V2 changes sign. In plain 
language, the former passes below and to the left of the latter for r = + and above and to the 
right for r — * +oo. By the Intermediate Value Theorem and continuous dependence, therefore, 
there exists at least one value ro for which they meet, i.e., there exists an undercompressive profile 
from the rest point associated with v% to the rest point associated with vi- 

Remark 4.6. The above, formal arguments, may be made rigorous as done in |FS] for the 
general nonisentropic case. They give slightly less information in the planar case (note that we 
lose the monotonicity /uniqueness of r* information obtained by phase plane analysis) but have the 
advantage of applying also to more general, nonplanar situations. 

4.3.2. The undercompressive bifurcation. There is an interesting bifurcation as r = /i/r decreases, 
between the situation of case (iii) in which there is a family of overcompressive connections between 
the rest points at v\ and t>4, bounded by Lax connections, and the situation of case (i), in which there 
exist no intermediate shock connections. This occurs at the point r = r* where an undercompressive 
connection appears. As illustrated in Figures [5] and [6j this occurs through squeezing of the infinite 
overcompressive family to a single under compressive-L ax profile pair, after which, as r is decreased 
past r*, the undercompressive connection breaks, leaving only the regular Lax connection and no 
intermediate profiles remain. 

Note that this occurs for the example in Figure [6] for value r* = 0.17, substantially less than the 



"physical" value predicted by (2.3) of r = /i/(2/u + r/) = 0.75. For the value r = .75 and 7 = 5/3, 7/5 



(monatomic or diatomic gas), we find numerically that undercompressive shocks do not occur. 

4.3.3. Composite-wave limits. In the limit as r — ► r+, the intermediate Lax shock connecting 
the rest points at 1)3 and v\ approaches a "doubly composite wave" formed by an approximate 
superposition of the limiting undercompressive profile between the rest points at V3 and V2 and the 
Lax profile between rest points at V2 and v\ at value r = r*, separated by an interval of length 
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going to infinity as r — > r+ on which the solution is approximately equal to the value of the saddle- 
type rest point at t>2 to which it passes nearby. Likewise, as r — > r+, the family of intermediate 
overcompressive profiles connecting the rest points at V4 and t>i approaches a triply composite wave 
consisting of the approximate superposition of the limiting Lax profile between rest points at 
and t>3, undercompressive profile between rest points at v% and V2, and Lax profile between rest 



points at V2 and v\ at value r 



separated by intervals of length going to infinity on which the 



solution stays near the saddle-type rest points at ^3 and V2- 

In either case, because the resulting profiles require larger and larger intervals in x to converge 
to limits U±, both the profiles and their associated Evans functions are numerically impractical 
to compute, requiring larger and larger computational domains, and must be handled separately 



taking into account the underlying limiting structure. We discuss this issue in Section 6.3 
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(a) 



(b) 



FIGURE 4. Phase portraits in singular limits. Figure (a) /i = 1, r = .1, Figure (b) 
H = .005, t = 1; parameters 7 = 5/3, v + = 0.1, / = 0.7, B2+ = 0.7. 
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(a) (6) 

Figure 5. Transition to nonexistence: keeping r = 2/x + 77 = 1 and letting /j — ► 0, 
we find that the overcompressive family is squeezed more and more but still connects 
until somewhere between fi = 0.185 (Figure (a)) and \i = 0.16 (Figure (b)). At that 
point the flows switch sides so that neither overcompressive nor undercompressive 
connections exist; parameters 7 = 5/3, v+ = 0.1, I = 0.7, B2+ = 0.7. 
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Figure 6. Undercompressive connection: at the transition point ~ \i = 0.17, an 
undercompressive profile appears; parameters 7 = 5/3, v + = 0.1, / = 0.8, B + = .7. 



5. Evans function formulation 



5.1. Two-dimensional MHD. In the two-dimensional case U3 = B% = 0, (2.16) becomes 

vt + v x - ui x = 
uu + u\ x + (aw -7 + Bl/(2no))x = t(ui x /v) x 

vB 2 ) t + (vB 2 ) x - Iu 2x = (crHoy 1 (B 2x /v) x , 



(5.1) 



where r = 2/i, + r/. Linearizing about the profile solution (v, ui,u 2 , B 2 ) we have 

vt + v x - ui x = 
uu + ui x + (-erytT 7-1 ?; + (1/ 'fio)(B 2 B 2 )) x = t(u\ x /v - u\ x v/v 2 ) x 

(5.2) 

u 2t + u 2x - (I/fi )B 2x = fi(u 2x /v-u 2x v/v 2 ) x 



at 



+ a x - Iu 2x = (0710) l (B 2x /v - B 2x v/v 2 ) x , 



where a = i)B 2 + vB 2 , so that B 2 = {a — B 2 v)/v. Substituting for B 2 we obtain the eigenvalue 
problem 



(5.3) 



where 



(5.4) 



Xv + v — u'i = 
Ami + u[ - (h(v)v/v l+1 )' = -{B 2 {a - B 2 v)/{fi v))' + t(u' x /v)' 



Xu 2 + u 2 — (I/no)(a/v — B 2 v/v)' = fi(u' 2 /v — u' 2 v/v 2 )' 



\a + a — Iu 2 = (crfio) (v (a/v — B 2 v/v)' — B' 2 v jv )', 



h(v) = —v' 1+1 (tu 1 /v 2 — a'yv 7 1 ) 



■P +1 (tv'/v 2 -ajv^- 1 ) 
-v~« + \v- 2 {v(v - 1) + av 1 - 1 -av+ {2^v)-\{B 2 ^ + Iu 2 ) 2 - v 2 B 2 _)) - ajv^- 1 ). 
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(a) 




(v,w)=(0.27,1 .24) 



(v ,w )=(0. 1,0.86) 



(v,w)=(0.85,-0.16) 
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(c) 

FIGURE 7. Typical phase portrait for MHD with two variables and infinite electric 
resistivity (a = oo). Parameter values are 7 = 2, v + = 0.1, I = 0.7, B 2+ = 0.7, 



and no = 1. In Figure (a) we graph f(v) given by (3.15). The Rankine-Hugoniot 
solutions corresponding to the roots in Figure (a) are given in Figure (b). In Figure 



(c) we plot level sets of <p(v, w) given by (2.42) and in Figure (d) we draw the phase 
portrait. 



(b) 



(d) 



(5.5) 



We let u(x) = f^ oo ui(z)dz, w = f^ oQ U2(z)dz, V = f_ v(z)dz, and a = J^ oQ a(z)dz to 
transform to integrated coordinates. Substituting we have 

XV + V" - u" = 

\v! + u" - (h(v)V'/^ +1 )' = -{B 2 (a' - B 2 V')/(fx v)y + r («"/«)' 

Xw' + w" — (1/ no)(a' /v — B2V fv)' = h{ w " /v — u^V/v 2 ) 1 

Xa' + a" - Iw" = ( ^ J M ) _1 (* _1 ((« , " B 2 V')/v)' - & 2 V/v 2 )'. 

Integrating from — cxd to x we obtain 

XV + V' - v! = 

Xu + u' - h(v)V'/v l+1 = -B 2 {a' - B 2 V')/{fi v) + t{u"/v) 

Xw + w' — (I / no)(a' /v — B 2 V /v) = fi{w" /v — u 2 V /v 2 ) 

Xa + a' - Iw 1 = (aiMi)~ x (v-\al fv - B 2 V' fv)' - B 2 V'/v 2 ). 
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(5.6) 




(a) 
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(c) 



(v ,w )=(0.1,0.21) 



-0.2 0.2 0.4 0.6 0.8 1 1.2 1.4 




(d) 



Figure 8. Typical phase portrait for MHD with two variables and infinite electric 
resistivity (a = oo). Parameter values are 7 = 5/3, v + = 0.1, I = 0.3, B2+ = 0.7, 



and no = 1. In Figure (a) we graph f(v) given by (3.15). The Rankine-Hugoniot 
solutions corresponding to the roots in Figure (a) are given in Figure (b). In Figure 



(c) we plot level sets of <p(v, w) given by (2.42) and in Figure (d) we draw the phase 
portrait. 



We use the coordinates (u, V, V, w, H w \ a i a> ' I '{° f° r the Evans function formulation. Solving 
for the desired variables, and using u" = XV + V", we have, finally, 



(5.7) 




cr/iov (X\iqV crfioT a/iQTV cr^ov 2 a^oTv' 1+1 



V. 
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(a) 

Figure 9 . Phase portrait in the planar case for MHD with two variables and infinite 
electric resistivity (a = oo). Parameter values are 7 = 5/3, v+ = 0.1, I = 0.7, 



(b) 



B 2+ = 0, and /j,Q = 1. In Figure (a) we plot level sets of <p(v, w) given by (2.42 ) and 
in Figure (b) we draw the phase portrait. 



This may be written as a first-order system W' = A{x, X)W from which the Evans function may 



be computed as described in Section 2.5 where 



(5.8) 
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(5.9) 
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2i?2 XB2 B 2 



R3 
°1 



B 2 v' B 2 h(v) 



. <T/^0^ C^O^ 07^)7" (THqTV CTfiQV 2 afJ,QTV^ +1 

W = (u, v, v', w, /xit/, a, a' / (a/j,Qv)) T , and /(■&) = r _1 (i) — f)™ 7 /i({>)). 



5.1.1. T/te case a = 00. In the case K3 = S3 = 0, with a = 00, the integrated eigenvalue equation 



(5.6) becomes 



(5.10) 



XV + V - v! = 
Xu + v! - h{v)V'/v^ +1 = -B 2 {a - B 2 V')/(fxov) + t(u"/v) 
Xw + w — (I / [io)(a /i) — B 2 V' /v) = fi(w" jv — u' 2 V'/v 2 ) 
Xa + a — Iw = 0. 



We use the coordinates (it, v, v', w, fiw', a) T . Solving for the desired variables, and using u" 
XV + V", K = I 2 / /i , we have 
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a' = Jw — Act. 

This may be written as a first-order system W' = A(x, X)W, where 



(5.12) 
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W = («,V,F',w,^«/ ) a)-', and f(v) = t-\v - v^h(v)). 

5.2. Three-dimensional stability. Finally, we consider the question of transverse stability, or 
stability with respect to three-dimensional perturbations, of a two-dimensional profile u 3 = B 3 = 0. 

Carrying third components through the computations of Section 5.1 we obtain the two-dimensional 
integrated eigenvalue equations (5.6) augmented with the additional equations 

(5.13) 



Xw 3 + w' 3 - (I/(j, )(a' 3 /v - B 3 V'/v) = n(w'l/v - u 3 V'/v 2 ), 

Xa 3 + a' 3 - Iw' 3 = (a/xor^Va/S " B 3 V'/v)' - B' 3 V'/v 2 ). 



In the case u 3 = B 3 = 0, these decouple from the rest of the equations, reducing to a transverse 
system 

Xw 3 + w 3 - (I/nov)a 3 = nw'l/v, 

Xa 3 + a 3 - Iw' 3 = {afiQv)~ 1 {a' 3 /v)\ 



(5.14) 



that may be studied separately. This is similar to the situation of the parallel case U2 = u 3 = B2 
B 3 studied in [FTllBiTZ] . 



5.2.1. The case a = 00. For a = 00, the transverse equations (5.14) become 

Iw' 3 - Xa 3 , 



a 3 



(5.15) 



IX 



fj,w 3 = Xvw 3 + vw' 3 — Kw' 3 H 03. 

Mo 

This may be written as a first-order system W' = A(x, X)W, where 

(5.16) A(x,X) = 

and W = (w 3 /jlw' 3 a 3 ) T , and used to compute a transverse Evans function determining stability 
with respect to perturbations in components u 3 , B 3 . 
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5.3. Construction of the Evans function. As described in |MaZ3[ IZlj . the above procedure 



may be carried out for general hyperbolic-parabolic systems under the standard assumptions (2.45 ) 



(2.46), (2.47), and (2.48), to express the eigenvalue problem as a first-order system W' = A(x, X)W 
with exponentially converging coefficient 

(5.17) \A(x,X)- A±(\)\ < Ce-^l, x ^ 0, 

where the constant C > is uniformly bounded on bounded domains in A. 

Moreover, defining 0, := {5RA > 0}, we have under the same standard hypotheses the general 
fact [MaZS^ that, on £1 \ {0}, the system satisfies the consistent splitting hypothesis of |AG Jj : the 
limiting coefficient matrices A±(\) have no center subspaces, and the dimensions of their stable 
and unstable subspaces agree and (by homotopy, using absence of center subspace) are constant 
throughout f2 \ {0}. Further, the associated eigenprojections, analytic on \ {0} by spectral 
separation (again, absence of center subspace), extend analytically to A = 0, so are analytic on all 
of the simply connected set O. By a standard construction of Kato [ Kato] , there exist analytically 
chosen bases (R^ , . . . , R^)(X) and (R^ +l , ■ ■ ■ R~ + N)(X) of the unstable subspace of A_ and the 
stable subspace of A + , respectively. 

Appealing to the general construction of Appendix [D] we may thus define the Evans function as 

(5.18) D(X) := det(Wf, . . . , , W+ +1 , . . . , W+)\ x=0 , 

where, for A 6 f2\{0}, {W- 1- } and {l^j"} are analytically chosen bases of the manifolds of solutions 
of W = A(x, X)W decaying as x — > +oo and x — > — oo, respectively, with 

(5.19) Wf(\,x) ~ e A± M x Rf(X) as i^ioo. 

Evidently, (i) D is analytic for A £ O, and (ii) A £ fl \ {0} is an eigenvalue if and only if D{X) = 0; 
sec |MaZ3j IZlj for further discussion. Moreover (see Appendix [Pj) , D is continuous with respect 



to model parameters (a, J, K) or {y+, B + , I). The asymptotics (5.19) may be used as the basis for 



numerical approximation of D; see [BrTl lBr2l IBrZl IBDG1 IHuZ2l IZ31IZ4] . 

6. Analytical stability results 

We begin by recording some analytical stability analyses in special cases, in particular certain as- 
ymptotic limits (small-amplitude, composite-wave, and high-frequency limits) that are numerically 
difficult. 

6.1. Small-amplitude stability. Consider first the small-amplitude limit, without loss of gener- 
ality (after rescaling) V- = 1 and v+ — * 1. 

Proposition 6.1. For J, K such that is a simple, genuinely nonlinear characteristic speed of the 
inviscid system for v+ = 1, there is a unique Lax-type profile connecting the rest points associated 
with v = 1 and v = v+ for v + sufficiently close to 1, and this profile is Evans, hence linearly and 
nonlinearly, stable (both with respect to coplanar and transverse perturbations). 



Proof. The existence result follows by a more general result of Pego |Pe] obtained by center manifold 
reduction. The stability result follows by a more general stability result of |HuZl] . □ 

6.2. Transverse stability of monotone profiles. Similarly as observed in the parallel case in 
[BHZj, in the case a = oo, transverse stability holds automatically for profiles that are monotone 
decreasing in v. Thus, in our numerical stability study, it is necessary to test transverse stability 
only for nonmonotone profiles. 

Proposition 6.2 ([BHZJ). For a = oo, monotone- density profiles, v x < 0, are Evans stable with 
respect to transverse perturbations: that is, they are three- dimensionally Evans stable if and only if 
they are two-dimensionally Evans stable. 
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Proof. Dropping subscripts, we may rewrite (5.15) in symmetric form as 



fiQvXw + fiovw — la = hhqw" , 
Xa + a — lw = I). 

Taking the real part of the complex L 2 -inner product of w against the first equation and a against 
the second equation and summing gives 



»A( J (vfi \w\ 2 + \a\ 2 ) = - J M/^ok'l 2 + J 



v x \w\ 2 < 0, 



a contradiction for 3?A > and w not identically zero. If w = on the other hand, we have a 
constant-coefficient equation for a, which is therefore Evans stable. □ 

Remark 6.3. In the above proof, we are using implicitly the fact that vanishing of the Evans 
function on 3?A > 0, A ^ 0, away from essential spectrum of L, implies existence of an eigenfunction 
decaying as x — ► ±oo |GJH IGJ2] , while vanishing at the point A = embedded in the essential 
spectrum implies existence of an L°° eigenfunction jZHt [MaZ3 . 

Remark 6.4. In the one-dimensional, parallel, case P>2 = -B3 = 0, U2 = u% = 0, the same argument 
yields stability if and only if the corresponding gas-dynamical shock is stable [BHZJ. (Recall that 
in this case MHD profiles reduce to gas dynamical profiles in (v,u\) |FT1 IBHZ] .) 

6.3. Stability of composite waves. We consider next the numerically difficult situation as de- 



scribed in Section 4.3.3 of a family of profiles passing closer and closer to one or more intermediate 
rest points, i.e., composite wave consisting of the approximate superposition of two or more com- 
ponent profiles separated by a distance going to infinity. This requires a computational domain 
[— L_,L+] of size going to infinity, hence arbitrarily large computational effort to resolve directly. 
However, it may be treated in straightforward fashion by a singular perturbation analysis taking 
account of the limiting structure. 



6.3.1. Double Lax configuration. Consider first the simplest case noted already in [Brl, Br2 of a 
family of overcompressive profiles in a four rest point configuration, bounded by two Lax/intermediate 
Lax profile pairs. Considering a family of overcompressive profiles U e connecting U- and U+ and 
passing closer and closer to an intermediate saddle £/*, parametrized by the distance e of the profile 
from U* , we find that the profiles approach composite waves consisting of the approximate super- 
position of the bounding Lax profiles U l and U 2 connecting U- to U* and E7* to U + , separated by 
a distance d(e) going to infinity as e — ► 0. 

Proposition 6.5. For the double-Lax configuration described, stability of U 1 and U 2 implies sta- 
bility of U e for e > sufficiently small. 

Proof. By a standard multi-wave argument, as described in the conservation law setting in |Z7| and 
(in slightly different periodic context) [OZ], the spectrum of any such composite wave U e approaches 
the direct sum of the spectra of its component waves U 1 and U 2 as e — > 0. More precisely, the 
Evans function D e associated with U e approaches a nonvanishing analytic multiple of the product 
of the Evans functions D 1 and D 2 associated with U 1 and U 2 , and so the zeros of D e approach the 
union of the zeros of D 1 and D 2 ; see |Z7] for details. In the case that U 1 and U 2 are stable Lax 
waves, D 1 and D 2 are nonvanishing on 5RA > 0, and so the the union of their zeros is empty. It 
follows that D e for e sufficiently small is nonvanishing on 3ft A > 0, giving the result. □ 
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6.3.2. Undercompressive configurations. Next, consider the more complicated examples of Section 
4.3.3 of doubly composite Lax profiles U e composed of the approximate superposition of a Lax 



profile U 1 and an undercompressive profile t/ 2 , and of triply composite overcompressive profiles 
U €1 ' e2 composed of the approximate superposition of a Lax profile U , an undercompressive profile 
U 2 , and a Lax profile [7 3 , where the parameter e (resp. e = (ei, €2)) indexes distance of the profile 
from the intermediate rest point (resp. points). 

Proposition 6.6. For either the Lax-undercompressive or Lax-undercompressive-Lax configura- 
tions described, stability of the component waves U 3 implies that D e has at most one unstable root. 



Proof. This follows by the observation, as in the proof of Proposition 6.5 that the zeros of D e 
approach the union of the zeros of D 3 as e — > 0, together with the fact that for stable Lax waves 
the associated Evans function has no zeros on 5RA > 0, while for stable undercompressive waves, 



the associated Evans function has a single zero at A = (see Proposition 2.8). From this we may 



conclude that D e has at most one zero on 3ftA > 0, giving the result. □ 

With the reduction to at most a single root, stability could in principle be decided as in 
[CHNZl IZ7| by examination of the mod two stability index of |GZ| IMaZ3] , a product T = 7 A 
of a transversality coefficient 7 for the traveling wave connection and a low-frequency stability de- 
terminant A, both real-valued, whose sign determines the parity of the number of unstable roots. 
The boundary case T = corresponds to instability through an extra root at A = [ZH, MaZ3 ; 
hence, Evans stability implies nonvanishing of T, 7, and A. 

The transversality coefficient 7 is a Wronskian of the linearized traveling- wave ODE measuring 
transversality of the intersection of the unstable manifold at [/_ with the stable manifold at [/+ 
of the traveling-wave ODE, with 7 / corresponding to transversality. From the composite wave 
structure, we may deduce that transversality of the component waves (a consequence of Evans 
stability, as noted above) implies transversality of U e for e > sufficiently small, or nonvanishing 
of 7; with further effort, the sign of 7 may be deduced as well. 

The low-frequency stability determinant A is for Lax shocks equal to the Lopatinski determinant 
determining inviscid stability; see the discussions of [ZS, Zl . In particular, it is independent of 
the nature of the viscous regularization, and readily computable. For overcompressive shocks, it 
involves also certain variations associated with the linearized traveling-wave ODE, as described 
in |Z3J IZlj . which though more complicated can also be computed in the limit e —* 0, deciding 
stability. 

In the Lax-undercompressive case, we can avoid such computations by the following observation. 

Corollary 6.7. For the Lax-undercompressive configurations, suppose that the limiting endstates 
U± = lim e ^o U± have a stable connecting viscous profile for some choice of viscosity ratios r = /i/r. 
Then, stability of the component waves U 3 , together with stability (resp. instability) of some 
composite wave U e ° for €q > sufficiently small implies stability (resp. instability) of all composite 
waves U e for e > sufficiently small. 

Proof. Since the composite wave is of Lax type, A is independent of the viscous regularization. It 
follows that A / for e > sufficiently small if there is an Evans stable profile for some choice 
of viscous regularization connecting the limiting endstates U±. (Recall from the discsussion above 
that Evans stability implies A/0 Ma/3llZT|.) Since 7 7^ for e > sufficiently small, as observed 
previously, we thus have that r / for e > sufficiently small, and thus T is of fixed sign. It 
follows that either all profiles U e are stable for e > sufficiently small, or no profiles U e are stable 
for e > sufficiently small, yielding the result. □ 



Using Corollary |6.7| we may conclude by numerical evaluations of component waves and a sample 
of composite waves with e small but nonzero the stability of composite Lax-undercompressive waves 
in the numerically inaccessible e —* limit. 
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Remark 6.8. Supposing that both standard Lax and composite Lax waves composed of Lax- 
undercompressive waves have been determined to be stable, and viewing the Lax-undercompressive- 
Lax composites as the composition of Lax=Lax-undercompressive and Lax waves, we obtain 
the partial result that triply composite waves U £1 ' e2 are stable for ei > sufficiently small and 
< €2 < E(ei), where E > depends on e±. However, to obtain a full result, it appears that one 
must carry out the more complicated computations described in the introductory discussion above, 
and so we do not complete this case. 

6.4. The large-amplitude limit. We now consider behavior as shocks of different types approach 
their maximal amplitudes. As computed in Appendix [C| for four rest point configurations, taking 
without loss of generality v\ < Vi < K < < = 1, the minimal value of v\ is and the 
maximum value of vi is 



2 V 4 



7 / 7 2 

v( J, K) = K + - - \ — + J(l -K)<K. 



Likewise, the minimum value of V3 is 



v( J, K) = K + 1 + yj^ + J(l - K) > K. 

Thus, for fixed J > 0, K > 0, the maximum-amplitude Lax 1-shock connects the rest points 
associated with v 4 = 1 and V3 = v > K, and the maximum-amplitude intermediate Lax 1-shock the 
rest points associated with V4 = 1 and t>2 = v < K. The maximum-amplitude Lax 2-shock connects 
the rest points associated with v 2 = v and v% — * 0, and the maximum- amplitude intermediate Lax 
2-shock the rest points associated with V3 = v and v\ — * 0. The maximum-amplitude (intermediate) 
overcompressive shock connects the rest points associated with V4 = 1 and v\ — > 0. For each of 
these limits, also a — > 0. 

(Here and below, we refer to two-dimensional shock types.) 

Proposition 6.9. For fixed J, K, the Evans function associated with Lax 1-shocks or intermediate 
Lax 1-shocks converges in the large- amplitude limit, uniformly on compact subsets of > 0, to 
the Evans function associated with the zero-pressure limit a = 0. 

Proof. An immediate consequence of the general property of continuous dependence on parameters 
of the Evans function, so long as the profile remains noncharacteristic and v remains bounded from 
the value v = at which the pressure function becomes singular. Noting that a = is bounded from 
the values a*(J,K) > and A(J, K) > at which profiles become characteristic (see Appendix 



C.3), and that v, v ^ 0, we obtain the result. □ 



The important implication of Proposition 6.9 is that stability of 1-shocks may be assessed nu- 



merically by computations on a finite mesh, even in the large-amplitude limit. 

Conjecture. We conjecture that, similarly, the Evans function associated with Lax 2-shocks 
or overcompressive shocks converge in the large- amplitude limit v+,a — > to an Evans function 
associated with the zero-pressure limit a = 0. 

Motivation. In the parallel case J = 0, this was shown by a delicate asymptotic ODE analysis in 
|HLZ| IBHZ] . Our numerics (Section [7]) indicate similar behavior in the general case; moreover, the 
limiting structure of the equations is quite similar, suggesting that the proof of |HLZ1 IBHZ] might 
extend with further care to nonzero values of J. 

6.4.1. Large- amplitude limit for transverse equations. As observed in [BHZ] . the coefficient matrix 



A(x, A) for the transverse eigenvalue system (5.16) is smooth (indeed, linear!) in the profile variable 



v, hence we obtain convergence in the large-amplitude limit of the transverse Evans function by the 
standard property of continuous dependence of the Evans function on parameters, Appendix [D] so 
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long as the profile v converges uniformly exponentially to its endstates, independent of a > 0, as it 
does in the regular limit arising for Lax 1-shocks, and appears numerically to do for Lax 2-shocks 
and overcompressive shocks as well. Our numerics (Section [7]) indeed suggest convergence. 

6.5. The high-frequency limit. Finally, we recall the following high-frequency asymptotics es- 
tablished in |HLyZl|, which we will use in our numerical studies to truncate the computational 
domain in A. 



Proposition 6.10 (|HLyZl ). Let D be the (integrated) Evans function associated with a non- 



characteristic shock profile of (2.1) (with either a = oo or a finite). Then, for some constants C, 
a, 

(6.2) lim D(X)/e aXl/2 = C, uniformly on 5RA > 0. 

|A|— >oo 

In particular, D does not vanish for 5RA > and |A| sufficiently large. 



Proof. This was proved in |HLyZl for isentropic gas dynamics in Lagrangian coordinates by an 



argument using the tracking lemma of |MaZ3[ IPZ] , However, the same argument applies to general 



hyperbolic-parabolic systems satisfying the standard hypotheses (2.45), (2.46), (2.47), (2.48), with 
the additional property that convection in hyperbolic modes is at constant speed. In this case, 
hyperbolic modes are specific volume v and, when a = oo, magnetic field B, each of which in 
Lagrangian coordinates are convected with constant speed — s = 1. Thus, the hypotheses are 
satisfied, and the result follows. (In the general case, D(X) ~ Ce aX ' for some a, (3, C.) □ 

7. Numerical stability investigation 

In this section, we discuss our approach to Evans function computation, which is used to de- 
termine whether any unstable eigenvalues exist in our system. Our approach follows the polar- 
coordinate method developed in |HuZ2] ; see also |BHRZ| IHLZ| |HLyZl IBHZ] . Since the Evans 



function is analytic in the region of interest, we can numerically compute its winding number in the 
right-half plane around a large semicircle B(0, R) D {3ftA > 0} appropriately chosen, thus enclos- 
ing all possible unstable roots. This allows us to systematically locate roots (and hence unstable 
eigenvalues) within. As a result, spectral stability can be determined, and in the case of instability, 
one can produce bifurcation diagrams to illustrate and observe its onset. This approach was first 
used by Evans and Feroe |EF] and has been applied to various systems since; see for example 
jPWllA^lBr^lBTJG] . 

7.1. Approximation of the profile. Following [BHRZ, HLZ] . we approximate the traveling wave 
profile using one of MATLAB's boundary-value solvers bvp4c |SGT] . bvp5c |KL| , or bvp6c [HM , 
which are adaptive Lobatto quadrature schemes and can be interchanged for our purposes. These 
calculations are performed on a finite computational domain [-L— , L + ] with projective boundary 
conditions M±(U — U±) = 0. The values of approximate plus and minus spatial infinity L± are 
determined experimentally by the requirement that the absolute error \U(±L±) — U±\ be within a 
prescribed tolerance, say TOL = 10~ 3 . For rigorous error/convergence bounds for these algorithms, 
see, e.g., |BeTl lBe2] . 

7.2. Approximation of the Evans function. Throughout our numerical study, we use the polar- 
coordinate method described in |HuZ2 , which encodes W = rfi, where 

W = W\ A • • • A Wk 

is the exterior product encoding the minors of W±, . . . ,Wk, "angle" $7 = u% A • • • A Uf. is the 
exterior product of an orthonormal basis {ujj} of Span{Wi, . . . , Wk} evolving independently of r 
by some implementation (e.g., Drury's method) of continuous orthogonalization and "radius" r is 
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a complex scalar evolving by a scalar ODE slaved to Q, related to Abel's formula for evolution 
of a full Wronskian; see [? HuZ21 \Z3\ IZ4] for further details. The Evans function is then recovered 
through 

£>(A) = WT A W + \ x=0 = det(Wf, . . . , W7, v W N )\ X=0 . 



Here, W 1 * 1 are approximated at x = —L^,L + using asymptotics (5.19) for by 



W~(-i_) ~ e'^ L -{R{ A---AR k ) 

where {-^J} is an analytically chosen basis for the unstable subspace U(A-) of A_ and /i = 
Traced- 1 {/(A-)) an d then evolved using the polar coordinate ODE toward the value x = where 
the Evans function is evaluated. The requirements on approximate plus and minus spatial infinity 
L± needed for accuracy are in practice the same as the requirement already imposed in the ap- 
proximation of the profile that the absolute error \U(±L±) — U±\ be within prescribed tolerance 
TOL = 10~ 3 ; see [HLyZl Section 5.3.4] for a complete discussion. L± = 10 sufficed for most 
parameter values. 

7.2.1. Shooting and initialization. The ODE calculations for individual A are carried out using 
MATLAB's ode45 routine, which is the adaptive 4th-order Runge-Kutta-Fehlberg method (RKF45). 
This method is known to have excellent accuracy with automatic error control. Typical runs in- 
volved roughly 60 mesh points per side, with error tolerance set to AbsTol = le-8 and RelTol 
= le-6. To produce analytically varying Evans function output, the initializing bases {Rf} are 
chosen analytically using Kato's ODE; see |GZ[ IHuZ2[ IBrZ[ IBHZ| for further discussion. Numeri- 
cal integration of Kato's ODE is carried out using a simple second-order algorithm introduced in 
IZ4| , a generalization of the first-order algorithm of |BrZ| . 



7.2.2. Winding number computation. We compute the winding number of the integrated Evans 
function D around the around the semicircle 

S :=d(B{0,R)n{^\ > 0}) 

by varying values of A along 20 points of the contour S, with mesh size taken quadratic in modulus 
to concentrate sample points near the origin where angles change more quickly, and summing the 
resulting changes in arg(Z)(A)), using 3dog-D(A) = argD(A)(mod27r), available in MATLAB by 
direct function calls. As a check on winding number accuracy, we test a posteriori that the change 
in D for each step is less than 0.2, and add mesh points, as necessary, to achieve this. Recall, by 
Rouche's Theorem, that accuracy is preserved so long as relative variation of D along each mesh 
interval remains less than 1.0. In Table [T] we give as a triple the radius of the domain contour, the 
number of mesh points, and the relative error for change in argument of -D(A) between steps. 
Care must be taken to choose R sufficiently large to ensure any unstable eigenvalues lie inside 



the domain contour S. Recall, Proposition 6.10 that 



(7.1) lim = C uniformly on 5RA > 0, 

v ; |A|— >oo e a * 



where a and C are constants. The knowledge that limit (7.1) exists allows us to determine a, C 
by curve fitting of log-D(A) = logC + aA 1 / 2 with respect to z := A 1 / 2 , for |A| >> 1. When D is 
initialized in the standard way on the real axis, so that D{\) = D(A), a and C are necessarily 
real. We then determine the necessary size R of the radius by a convergence study, taking R to be 
a value for which the relative error between D(\) and Ce a ^ becomes less than .1 on the entire 
semicircle with 3ft A > 0, indicating sufficient convergence to ensure nonvanishing. (Relative error 
< 1 implies nonvanishing.) For many parameter combinations, R = 2 was sufficiently large, though 
some required a much larger radius. 
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Remark 7.1. Alternatively, we could use energy estimates or direct tracking bounds as in [HLZ 
and |HLyZl| , respectively, to eliminate the possibility of eigenvalues of sufficiently high frequency. 
However, we have found the convergence study to be much more efficient in practice; see [HLyZl| . 

7.3. Description of experiments: broad range. In our numerical study, we covered a broad 
intermediate parameter range to demonstrate stability of Lax and overcompressive profiles. To 
avoid redundancy, we discarded four rest-point configurations for which v = 1 was not the largest 
(v- value of a) rest point, since these can always be rescaled to an equivalent configuration for which 
v = 1 is largest, hence otherwise would be counted twice. The following parameter combinations 
were examined, when physically meaningful, for Evans stability: 

(j,v + ,I,B2+,iMi) € {7/5,5/3} 

x {0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, HT 1 , 10" 2 } 
x {0.2, 0.4, 0.6, 0.8, 1.2, 1.4, 1.6, 1.8, 2.0} 
x {0.2,0.4,0.6,0.8, 1.0, 1.2, 1.4, 1.6, 1.8,2.0} 
x{1.0}. 

For w + = 10 -2 above, the Mach number, as computed in appendix [bJ typically varies between 
20 and 40. For a little over 30 of the parameter combinations above for which / > 1, we took 
v + = 10~ 3 , 10 -4 , and 10 -5 attaining a Mach number of over 10,000 in some cases. All Evans 
function computations were consistent with stability. 

We also covered a broad intermediate range in terms of the parameters (K, J, a) . When physically 
relevant we examined the parameter combinations: 

(-y,K,J,v + ,no) e {7/5,5/3} 

x {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 1.05, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0} 
x {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0} 
x {0.1} 
x {1.0}. 

Finally, we examined the stability of the whole family of over compressive profiles for the relevant 
parameters belonging to 

(-y,K,J,a,no) £ {7/5,5/3} 

x {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95} 

x {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0} 

x {01,02,03,04,05} 
x {1.0}. 



v+ 


B + = 0.2 


B + = 0.8 


B + = 1.2 


B + = 1.6 


B + = 2 


0.1 


(2,20,1.5(-1) 


(4,20,1.8(-1) 


(8,64,9.1(-2)) 


(16,64,1.3(-1) 


(2,20,7.6(-2)) 


0.4 


(2,20,1.1(-1)) 


(2,20,3.3(-2)) 


(2,20,3.1(-2)) 


(2,20,3.1(-2)) 


(2,20,3.2(-2)) 


0.6 


(2,20,7.7(-2)) 


(2,20,1.6(-2)) 


(2,20,1.8(-2)) 


(2,20,1.7(-2)) 


(2,20,1.8(-2)) 


0.8 


(2,20,6.6(-2)) 


(2,20,2.0(-2)) 


(2,20,1.9(-2)) 


(2,20,1.9(-2)) 


(2,20,2.0(-2)) 



Table 1. Table demonstrating contour radius, number of mesh points, and relative 
error. Here / = 1.2 and 7 = 5/3. 
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where a\ = 10~ 3 and 05 is the largest value of a such that the system has 4 fixed points of the form 
(y,w) with v < 1. For each value a,- we examined the stability of 5 profiles chosen by requiring 
they pass through evenly spaced points along the line in the phase plane connecting the two rest 
points with intermediate v+ coordinates, thus insuring our profiles be representative of the family 
of over compressive traveling waves. In Figure 10 we plot in bold some profiles examined in our 
over compressive study. 




(o) (6) 

Figure 10. The bold curves in the phase portrait are the over compressive profiles 
for which the integrated Evans function -D(A) was computed. The parameter values 
for Figure (a) are 7 = 5/3, v + = 0.1, I = 0.7, B2+ = 0.7, and fiQ = 1. In Figure (b) 
we have 7 = 5/3, v + = 0.1, I = 0.6, -B2+ = 0.9, and //q = 1- 



In the whole investigation, each contour computed consisted of at least 40 points in A. In all 
cases, we found the system to be Evans stable. Typical output is given in Figure [TTJ We remark 
that the Evans function is symmetric under reflections along the real axis (conjugation). Hence, 
we only needed to compute along half of the contour (usually 20 points in the first quadrant) to 
produce our results. 




0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



(a) (6) 

Figure 11. Typical Evans function output. The parameter values for Figure (a) 
are 7 = 5/3, v + = 0.1, / = 0.7, B2+ = 0.7, and fio = 1- I n Figure (b) we have 
7 = 5/3, v + = 0.1, / = 1.4, B 2+ = 1.4, and ^ = 1- 
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7.4. Composite limit. As described in Section 6.3 as overcompressive shocks approach the limit 
of a composite wave formed by the approximate superposition of the bounding Lax 1 and 3-shocks, 
separated by larger and larger distance, the Evans function computation becomes prohibitively 
costly. However, the analytical result of Proposition |6.5| shows that we need not carry that out, 
since stability in the composite limit follows by stability of the component Lax waves, already 
tested. 



7.5. Large-amplitude limit. As shown analytically in Section 6.4 the Evans functions for Lax 
1-shocks and intermediate Lax 1-shocks converge in the large-amplitude limit a — > 0, both for 
coplanar and transverse perturbations, in this case the nonphysical boundary a — > need not be 
treated in any special way. 

We carried out numerical case studies suggesting that the Evans function converges in the large- 
amplitude limit also for the more singular cases of Lax 2-shocks, intermediate Lax 2-shocks, and 



intermediate overcompressive shocks, left unresolved in the analytical treatment of Section 6.4 
We conjecture that convergence holds also in these cases, as shown in the parallel case J = in 
[HLZllEHZ] . 

A case study of the Lax 2-shock and intermediate Lax 2-shock cases is displayed in Figure 
12 corresponding to a two-rest point configuration, with parameters K = 2, J = 1, 7 = 5/3, 
a = 10~ 3 , 10~ 4 , 10~ 5 , 10~ 6 , 1CP 7 , 10~ 8 . We found stability for all amplitudes in each of these cases. 
We only had to take contour radius R = 16 for a as small as a = 10 -8 , so these runs were not 
computationally expensive. The Mach number for a = 1CP 8 is ~ 10, 954. 




0.2 0.4 



(a) 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



(b) 



Figure 12. Large-amplitude limits, parameters K = 2, J = 1, 7 = 5/3. In Figure 
(a), we display the image of the semicircle under D for a Lax 2-shock in two-rest 
point configuration in the a — > limit, a = 10~ 3 , 10~ 4 , 10 -5 , 10 -6 , 10 -7 , 10~ 8 . where 
a = 10~ 8 corresponds to Mach number ~ 10, 954. Convergence of contours appears 
to occur at a ~ 10 -6 . or Mach number ~ 1,095. In Figure (b), for the same 
sequence of a-values, we display the images under the transverse Evans function, 
again suggestive of convergence. 

A case study of the Lax 2-shock, intermediate Lax 2-shock, and intermediate overcompressive 
cases is displayed in Figure [13] corresponding to a four-rest point configuration, with parameters 
K = 0.7, J = 0.5, and a = 10 _1 , 10" 2 , . . . , 10" fc , taking a as small as necessary to achieve con- 
vergence: for example, in the overcompressive case, a = 10 -7 , or Mach number ~ 3,817. In each 
case, convergence was achieved; likewise, we again found stability for all amplitudes. See Figure 
for the corresponding phase portrait with K = 0.7, J = 0.5, and a = 10 -8 ~ 0, approximating 
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the a — ► limit. Note that each of the Lax 2-shock, intermediate Lax 2-shock, and intermediate 
overcompressive shock profiles appear to lie on a straight line orbit. It would be interesting to check 
whether the a = traveling-wave ODE, a polynomial (cubic) vector field, indeed supports exact 
straight line connections. 




! 0.2 0.4 0.6 0.8 1 1.2 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

(c) (d) 

Figure 13. Large amplitude limits, parameters K = 0.7, J = 0.5, and a = 
10 _1 , 10~ 2 , . . . , 10~ fc , getting smaller as necessary to see what appear to be con- 
vergence to a limit, (a). Lax 2-shock, vi to V\. (b). Intermediate Lax 2-shock, v% 
to v\. (c). overcompressive 1-2 shock, v<± to v\. (d). Transverse Evans study for (c). 
In each case, we appear to obtain convergence at a = 10~ 7 , corresponding to Mach 
number ~ 3, 817. 



7.6. Three-dimensional stability. As discussed in section [6] transverse stability holds automat- 
ically in the case a = oo for profiles that are monotone decreasing in v, so for all the studies 
described previously, we examined stability of (5.15 ) only in the case of a nonmonotone profile. All 
computations were consistent with stability. 
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Figure 14. Phase portrait corresponding to Figure 13 parameter values K = 0.7, 
J = 0.5, and a = 10~ 8 ~ 0. 




-20 -15 -10 -5 




Figure 15. Typical transverse Evans function output, parameter values 7 = 5/3, 
I = 0.6, B + = 1.4, and /j,q = 1. In Figure (a) we display the nonmonotone profile. 
In Figure (b) we display the winding number computation. 

7.7. The undercompressive case. For our undercompressive study we, taking care to avoid 
repetitions, considered the parameter combinations 

^,v+,K,J) €{7/5,5/3} 

x {0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9} 
x {0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9} 
x {0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9} 

for which a four rest point configuration exists in phase space. We fixed 77 and let /1 be a free 
parameter in the boundary value problem allowing us to solve simultaneously the value of /i for 
which a traveling wave connects the saddle points and for the profile itself. We successfully found 
and examined over 250 undercompressive profiles for Evans stability. Since undercompressive pro- 
files are monotone implying stability in the transverse case, we only computed the Evans function 



associated with (5.12). Because the Evans function output for undercompressive waves has a zero 
at the origin, we used a small half circle of radius 10 -3 as part of the domain contour to skirt 
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around the origin. All winding number results were consistent with stability. Typical output is 
displayed in Figure 16 



0.2 - 
0.15 - 

0.1 - 
0.05 - 

- 
-0.05 - 

-0.1 - 
-0.15 - 

-0.2 - 




-0.05 
-0.1 

-0.15 
-0.2 

-0.25 




(a) 

FIGURE 16. Typical Evans function output, undercompressive case; parameter val- 
ues v+ = .3809; I = 0.8; B+ = 1.1692; 7 = 5/3; n = 0.2381; and 2/i + f] = 1. In 
Figure (a), our domain contour comes within 10 -4 of the origin, with a gap between 
mesh points ilO and — il0~ 4 , and our range contour comes within 10~ 6 of the ori- 
gin. In Figure (b), our domain contour comes within 10 -3 of the origin, then follows 
a small semicirle around it, the image of which may be seen at the far righthand 
side of the figure. (For undercompressive shocks, the integrated Evans function has 



a zero at the origin |ZH| IMaZ3j : see Proposition 2.8 ) 



8. Discussion and open problems 

In this paper, we have carried out by a combination of asymptotic ODE analysis and numerical 
Evans function computations a global existence/stability study for viscous shock profiles of two- 
dimensional isentropic magnetohydrodynamics with infinite electrical resistivity. For a monatomic 
7-gas equation of state, and standard viscosity ration ij = —2/^/3, we find that Lax and over- 
compressive profiles appear but undercompressive profiles do not. A systematic numerical Evans 
function investigation indicates that all profiles are nonlinearly stable both with respect to two- 
dimensional and three-dimensional perturbations. For different viscosity ratios, undercompressive 
shocks can appear, and these appear also to be stable with respect to two- and three-dimensional 
perturbations. 

Our stability analysis generalizes previous viscous studies of the viscous stability problem in 
\FT\ IBHZ] for the parallel case. See also the investigations of stability in the small-magnetic 
field limit in (MeZ, GMWZ2 . For analyses of the related inviscid stability problem, see, e.g., 
jTJ \BT\ IMeZ] and references therein. 

Much of the analysis carries over to the full three-dimensional case; in particular, the Rankine- 
Hugoniot analysis is completely general. It would be very interesting to carry out a systematic 
analysis in three dimensions, following the approach laid out here. Genuinely three-dimensional 
profiles, having a richer structure and more degrees of freedom, would appear to be a good place to 
look for possible instability or bifurcation. (Here, three-dimensional refers as in the present paper 
to the dependent variables and not the independent variable x.) As noted in [TZ , instability, 
by stability index considerations, would for an ideal gas equation of state imply the interesting 
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phenomenon of Hopf bifurcation to time-periodic, or "galloping" behavior at the transition to 
instability. 

Likewise, it would be interesting to carry out the full, nonisentropic case, building on the 
Rankine-Hugoniot study of [FRlJ. This should in principle be straightforward using the meth- 
ods developed here, but more computationally intensive. We suspect that, as in the gas-dynamical 
case |HLyZl , the large-amplitude limit may in fact be straightforward in this case, with w + bounded 
from the value zero at which the pressure function has a singularity. Further interesting general- 
izations would be to consider a Van der Waal or other "real gas" equation of state, or to include 
third-order dispersive effects modeling a "Hall effect" as in [DR . 

Interesting boundary cases left open in the present analysis are the K — > oo singular perturbation 



problem discussed in Section C.2 and the large-amplitude limit v+ — > for Lax 2-shocks and 
over compressives. We conjecture that in each of these cases, the Evans function converges in the 
limit to the Evans function for the formal pressureless gas limit a = 0. (As shown in Appendix [C] 
all of these limits coincide with a — ► 0.) Our numerics are consistent with this conjecture, while 
the related analysis of |HLZ| gives an idea how to prove it. At the physical, modeling, level, an 
important problem is to determine physically interesting values of J, K, a, and the viscosity ratio 
r = /i/(2/x + rj), which strongly affects solution structure as we have seen. 

The absence of instabilities in our experiments suggests perhaps the more general question 
whether shock profiles for systems possessing a convex entropy are always stable. We do not at the 
moment see why this should be so, and suspect that an ideal gas equation of state is perhaps too 
simple an example on which to base conclusions. However, to verify or produce a counterexample 
to this conjecture takes on a larger importance in light of the growing body of stable examples. 

Another interesting direction for further investigation would be a corresponding comprehensive 
study of multi-dimensional stability, as carried out for gas-dynamical shocks in |HLyZ2 . (Here, 
multi-dimensional refers to the independent variable x.) As pointed out in |FT| . instability results 
of |BXl [T] for the corresponding inviscid problem imply that parallel shock layers become multi- 
dimensionally unstable for large enough magnetic field, by the general result |ZS| \Z1\ IZ2| that 
inviscid stability is necessary for viscous stability, so that in multi-dimensions instability definitely 
occurs. The question in this case is whether viscous effects can hasten the onset of instability, that 
is, whether viscous instability can occur in the presence of inviscid stability. 

Appendix A. Signature of V 2 <f> 



Taking the Hessian of the relative entropy <p defined in (2.40), using (2.39), we readily obtain 



(A- 1 ) V^ )B/Mo) <£=| 1 -I 



p'iv) + 1 Ho{B/^) 

1 -I 

Ho(B/hq) -I vfi 



yielding after a brief computation 



det ^l, w , B /^ = Mp'(v) + l)(v — K) - \B\ 2 



At rest points of the traveling wave equation, we have B = J((l — K) /{v — K ) , by (3.9), yielding 
detV 2 „ 

w^/no)^ - IM)(v — K)f'(v), where / is as defined in ( |3.13| ). This gives an explicit connection 



between the signature of V 2 <^ and the sign of the derivative of the reduced Rankine-Hugoniot 
function /. 

Appendix B. Computing the Mach number 



The Mach number of a Lax shock is defined as 

u- — a 
M = , 

c_ 
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where u_ is the downwind velocity, a is the shock speed, and c_ is the downwind sound speed (in 
the characteristic family of the shock), all in Eulerian coordinates. (Here, the downwind side is at 
— oo, since we consider a left-moving shock.) By considering the conservation of mass equation, we 
have pt + (pu) x = 0, where p = 1/v is density. Hence, the jump condition is given by a[p] = [pu], 



which implies, in the original scaling for (2.1), that 

U + V- — U-V+ 



a 

V- — v + 



Hence, 



u- — a v-{u- — u+) v-[u] V- 
C_ c_(t>_ — v + ) C-[v] C- 



Noting that < v+ < u_ = 1, we simplify to get 

c_ 

where c_ is the sound speed in Eulerian coordinates at v— = 1, u__ =0, w_ = 0. 
Here, for Lax 1-shocks (resp. 2-shocks), 

<?- = \ ([eg + (Pto)-\B1 + I 2 )] ± ^[cl + { P p,)-\Bl + P)Y - ^(ppo)- 1 ! 2 



1 

2 

by (2.8), where c s = dp/dp\ p= i = y/dq denotes sound speed \K \ IMeZj . In the parallel case J = 0, 



K [ja + 2J + K]± ^[ja + 2J + K] 2 - 4jaK) 



this gives c_ = ^/ja, or M = —/=%■> f° r Lax 2-shocks or Lax 1-shocks with 7a > K, in agreement 



with the standard Mach number for gas-dynamical shocks. However, for Lax 1-shocks with 7a < K, 
it gives the anomalous value M = 1/K. (Recall that MHD profiles in the parallel case reduce to 
gas-dynamical profiles.) 



It is readily verified that the Mach number is invariant under the rescaling (2.15), hence gives a 
useful measure of shock strength in the original unrescaled coordinates. However, the example of 
the parallel case shows that it can give anomalous values for other than the simple gas-dynamical 
case. Note also that this measure of shock strength involves only the left state at — 00, so does not 
distinguish between intermediate vs. regular types of shocks. We therefore use the Mach number 
only to give a rough idea of the strength of shocks considered in our studies, and not as a systematic 
measure of strength across all parameters. 

Appendix C. Limiting cases for Lax and overcompressive shocks 



C.l. a — > 0. One may ask the questions of the limits of positivity of a given by (3.10). One has 
^r^^ l_-v+ ( 1 + V+-2K 

(C - X) a = ^iV- J {v + -K? 

The two roots of the term between brackets are real if and only if i > 1 — K, the bracket being 
always positive when J < 4(1 — K). In the former case, there are two roots v < v, with a < for 

v < v + < v. Moreover, v(J, K)=K + ±- ^ + J(1-K), so v(J, K) > for K 2 - J{l-2K) > 0. 

If K > g, this is always true. If K < i, it is true for J < jj^U , with v = precisely on the limiting 

curve J = a bove which there are no rest points with values v+ < K. 

Returning to the discussion, we have thus, for the two roots between and K, that the limiting 
values of v+ for which a — > are v+ — *■ 0, for which the factor goes to zero while the factor 

V, — 1 

in brackets remains bounded, and — » v > 0, with the two limits coinciding precisely in the 
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case J = i^njc This means, for fixed (J,K), that v+ — > in the large-amplitude limit for Lax 
2-shocks or overcompressives (the ones involving the rest point v\ with smallest v-value, and for 
which v\ < K), except on the measure zero boundary J = j^k- On the other hand, noting that 
K (since the term in brackets is negative there) lies always between v and v, values v+ for Lax 
2-shocks and overcompressives remain bounded away from the value v+ = K at which a becomes 
singular, so long as J 7^ 0. 

That is, for Lax 2-shocks and overcompressives, the boundary of existence for fixed J, K is 
marked by the singular limit v + — > on one side, and v + — ► v( J, K) > on the other side. In the 
measure-zero case J = > ^ ne ^ wo nm its agree. For Lax 1-shocks on the other hand, K < v+, 
hence v+ — ► only in the limit as J, K — > 0, and so this singular limit does not arise for J, K 
fixed. Neither does v + — > K, unless J = 0. 

C.2. a — ► 00? For K < 1, the upper limit for a is the characteristic boundary a < A(J, K) described 



in Proposition 3.11 For K > 1, may take without loss of generality a < a m i n (J, K), where a m i n is 
defined as the minimum value of a at which rest points v > K appear, since in the latter case one 
can then rescale to the case K < 1 already treated. Thus, for J > and < K ^ 1 bounded, a 
may be taken always finite. 

In the case K < 1, we've already seen in Remark 3.10| that J is finite for fixed K > bounded 



from 1. In the case K > 1, by Remark 3.12 we have J < 4(K — 1) to begin with, once we eliminate 



four rest point configurations (as we may do by rescaling so that K < 1). So, for fixed K > 1, we 



get bounded J, K, hence bounded a by (3.10) and < v+ < 1. Combining these observations, we 
find for any bounded K that is also bounded from 1 that J and a may be taken bounded as well. 
Thus, we need only consider finite parameter values (a, J) for K bounded and bounded from 1. 
The sole remaining case is that of a Lax 2-shock, v + < 1 < K, with K going to infinity and 



J < A(K — 1). Consulting again (C.l ), we see that in this case a~J/K<4&sK^ 00, so that a is 
again uniformly bounded. (On the other hand, this case can certainly occur for a sufficiently small.) 
The conclusion is that, without loss of generality (i.e., rescaling four rest point configurations to 
K < 1 whenever they occur), a may be taken uniformly bounded, independent of J, K, so that 
a — ► 00 does not occur. On the other hand, the case K — > 00 can occur, and even J — > 00, K — > 00 
simultaneously, with a remaining finite. 

In this latter case, the profile ODE, and the associated stability problem, should be treatable by 
a singular perturbation analysis, rescaling x, in which the pressure term disappears. However, we 
do not carry out this analysis here. 



C.3. Characteristic boundaries. Other important limits are the parameter values for which 
the shock becomes characteristic at U+ or [/_, since the rate of exponential convergence of the 
shock profile goes to zero as they are approached, so that the length of the computational domain 
[— L_,L+] needed for accurate numerical approximation goes to infinity. These are given by the 
surfaces a = a*(J, K) = ^rpf^v (corresponding to v + = 1) and a = A(J,K) (corresponding to 
vi = V2). The first is resolved by the analytical result of small-amplitude stability. The second 
requires a refined analysis outside of the scope of this paper, involving stability of characterstic 
shocks. For results in this direction, see |HoZ2l iHo] . 

Appendix D. The conjugation lemma 
Consider a general first-order system 
(D.l) W' = A p {x,X)W 

with asymptotic limits A^. as x — > ±00, where p £ M. m denote model parameters. 
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Lemma D.l ( [MeZHlFZ] ). Suppose for fixed 9 > and C > that 
(D.2) \A P - A p ± \(x,X) < Ce-*W 

for x ^ uniformly for (X,p) in a neighborhood of (Ao), Po and that A varies analytically in X and 
smoothly (resp. continuously) in p as a function into L°°(x). Then, there exist in a neighborhood 
of (Ao,Po) invertible linear transformations P p (x, A) = I + 0+(cc, A) and P_(a;, A) = / + Q p _(x, A) 
defined on x > and x < 0, respectively, analytic in A and smooth (resp. continuous) in p as 
functions into L°°[0,±oo), such that 

(D.3) \O p ±\ < C ie - §M forx^O, 

for any < 9 < 9, some C\ = C\(9,9) > 0, and the change of coordinates W =: P±Z reduces 
( |D.l ) to the constant- coefficient limiting systems 

(D.4) Z' = A P ± Z for x ^ 0. 

Proof. The conjugators are constructed by a fixed point argument |MeZlj as the solution of an 
integral equation corresponding to the homological equation 

(D.5) P' = A p P - A P ± P. 

The exponential decay (D.2) is needed to make the integral equation contractive with respect to 
L°°[Af, +oo) for M sufficiently large. Continuity of P± with respect to p (resp. analyticity with 
respect to A) then follow by continuous (resp. analytic) dependence on parameters of fixed point 
solutions. Here, we are using also the fact that ( |D.2 ) plus continuity of A p from p — > L°° together 
imply continuity of e^(A p — A p ±) from p into L°°[0, ±oo) for any < 9 < 9, in order to obtain the 
needed continuity from p — > L°° of the fixed point mapping. See also |PZ]. □ 

Definition D.2 (Abstract Evans function). Suppose that on the interior of a set £1 in X, p, the 

dimensions of the stable and unstable subspaces of A±(X) remain constant, and agree at ±oo ("con- 
sistent splitting" |AGJ] ). and that these subspaces have analytic bases extending continuously 
to boundary points of £1. Then, the Evans function is defined on £1 as 

(D.6) D p (X) := det(P+ J R+ . . . , P+P+ P~ R k+1 , p-R' N )\ x=0 , 

where P|! are as in Lemma 



D.l 



Evidently, Wf := P + P/, j = 1, • • • , k} and Wj := P~ Rj , j = k + 1, . . . , N} are bases for the 
manifolds of solutions W of W = A P W decaying as x — > +oo and x — > — oo, respectively, analytic 
in A and smooth in x, with Wf{x) ~ e A±x Rj± as x — > ±oo. (Here, we suppress p for notational 
convenience.) Thus, D p has the alternative representation 

D p (X) = det(iy+, . . . , W+, Wj +V W N ) X=0 



given in (2.56). 

Using the fact that P^ — > / as x — > ±oo, and the fact that modes W- are growing as x goes 
from ±oo to with undesired modes exponentially decaying, it is not difficult to see that the Evans 
function can be well-approximated by replacing Wj 1 with solutions w^ ,approx of (D.f ) with data 

w ±,appro^ ±L) = e A ± (X)(±L) R ± = p±(± L )-l W ±(± L ). 

(Compare solutions Z^ = P^Wp and z^' approx = p±w^> a PP rox Q f the constant-coefficient equa- 
tions Z' = A±Z.) This can be used as the basis for numerical approximation of the Evans function, 
as described, e.g., in [BED EEH EH iBDGl lHuZ2l IZ3llZ4] . 
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